Method and apparatus for wide field distortion-compensated imaging

ABSTRACT

An imaging system for measuring the field variance of distorted light waves collects a set of short exposure &#34;distorted&#34; images of an object, and applies a field variant data processing methodology in the digital domain, resulting in an image estimate which approaches the diffraction limited resolution of the underlying physical imaging system as if the distorting mechanism were not present. By explicitly quantifying and compensating for the field variance of the distorting media, arbitrarily wide fields can be imaged, well beyond the prior art limits imposed by isoplanatism. The preferred embodiment comprehensively eliminates the blurring effects of the atmosphere for ground based telescopes, removing a serious limitation that has plagued the use of telescopes since the time of Newton.

FIELD OF THE INVENTION

The present invention relates to the imaging of objects which have beendistorted or blurred, and more particularly relates to a method andapparatus for creating clear wide field images after firstcharacterizing and then compensating for the distortion. The distortionis characterized by being both spatially and time variant.

BACKGROUND AND SUMMARY OF THE INVENTION

The limitations on imaging system performance imposed by a turbulentmedia, most simply described as `blurring,` are well known, particularlyin applications using medium to large aperture telescopes in the openatmosphere. These limitations have not only led to a variety of systemsolutions that will be discussed as prior art, but have played a majorrole in the decision to launch space based telescopes and have led toserious postulations of lunar based observatories.

For a large aperture telescope--generally greater than a 10 centimeterdiameter for the visible light region--which is otherwise constructed toa precision commonly referred to as "near diffraction limited," theoverall ability to resolve objects obscured by a turbulent atmosphere islimited by the turbulence rather than by the instrument. For the visualband of light once more, it is quite typical for a 1 meter aperturetelescope to have ten times worse resolving power due to the turbulence,while a 10 meter aperture telescope can be 100 times or more worse thanits innate "diffraction limit." The exact numbers for any giventelescope on any given night are a function of many variables, but thisgeneral level of degradation is widely recognized. As importantly, thisatmospheric blurring directly leads to a loss in effective sensitivityof these large aperture imaging systems, which either renders dimobjects just too dim to be seen or forces greatly extended exposuretimes, ultimately limiting the number of objects that can be imagedduring a given length of usage time.

The prior art for addressing this problem and trying to alleviate it canbe generally categorized into the following well known areas: 1)Telescope Placement; 2) Adaptive Optics Systems; and 3) SpeckleInferometric Systems. It would not be unfair to say that the systemdisclosed herein is best summarized as a fundamental expansion to thethird category, though this is only in a most general sense.

Regarding the first category, that of simply finding locations ofminimal turbulence, the placement of telescopes at high elevations hasbeen known and practiced since Isaac Newton's area. This typically addssome expense, but more critically, it is well known that this only goesa small way toward reaching diffraction limited performance, andmoreover, imaging performance is still quite subject to the varyingatmospheric conditions of a given site. The space age has brought aboutthe "obvious" solution of launching telescopes above the atmosphere,through at considerable expense and complexity.

The second category of adaptive optics has been well known for decadesand has seen significant physical realizations over the last twodecades. A brief technical summary of such a system is that after atelescope primary mirror has collected the light waves emanating from agiven object, it splits the light wave into two "beams." One beam goesinto an instrument known as a wavefront sensor and the other beam entersan ordinary imaging system or camera. The wavefront sensor derivesinformation about the phase distortion (caused by the atmosphere) of theincoming light wave and in less than hundredths or even thousandths of asecond, sends control signals to mirror actuators which advance andretard primary beam mirror surfaces in cancelling opposition to theincoming phase distortion. There are two critical problems with thesesystems, however. First, they are expensive to build and expensive tomaintain. Second, they can only increase the resolving power within anextremely small angle about the nominal central (paraxial) ray of thetelescope, typically on the order of 2 to 10 arc seconds in the visibleband. This "tunnel vision" is technically referred to as the"isoplanatic patch." It arises due to the fact that the phase distortionof a three dimensional media such as the atmosphere changes rapidly as afunction of field angle.

The third category of speckle interferometry was begun and developedroughly contemporaneously with adaptive optics systems, i.e. over thelast two decades. The guiding principle of these systems is to take alarge number of short exposure images of an object, thereby "freezing"the phase distortion (in time) of a given single exposure, andultimately deriving the fourier magnitude and the fourier phase of theobject being imaged through mathematical operations on the accumulatedset of exposures. The state of the art of these techniques now includesthe addition of a wavefront sensor to an overall system, which providesadditional information to the computational operations, allowing for asmaller number of short exposures to be gathered to achieve a givenlevel of precision and also simplifying the overall process. Thoughthese speckle systems can be less expensive than adaptive opticalsystems, they too are limited by "tunnel vision" of the isoplanaticpatch. Moreover, since they do not account for the field variance of thedistortion, they suffer from higher noise and error levels, which inturn affect their effective sensitivity, accuracy, and throughput(number of images needed for a given level of accuracy).

The expense and limited performance of the collective prior art hasaccordingly limited its application and hampered its broader acceptance.A novel aspect of the system disclosed herein is that it will quantifyand compensate for the field variant complex phase distortion across anarbitrarily wide field of view, at a cost much less than a comparableadaptive optics system, and with a throughput and accuracy much betterthan speckle interferometric systems. In so doing, it can also provideadvantages over these prior art systems even in applications where theobject of interest is contained within the so-called isoplanatic patch;most notably in cost per image at a given accuracy (error and signal tonoise ratio).

Just as adaptive optics systems have recently employed "artificialbeacons" to assist in the imaging of very dim objects, so too can thesystems described herein utilize various forms of this concept asdescribed herein. Artificial beacons can be employed when the brightnessof an object under study is insufficient or inappropriate to providephotons to a wavefront sensor. The beacons are generally laser beamsdirected along a close line of sight to the object, generatingbackscatter photons which will undergo largely similar phase distortionsas the photons from the object under study, and thus they can be used todeduce the phase distortions applicable to the object photons.

The preferred embodiment of the system disclosed herein is a largetelescopic imaging system used in the open atmosphere. This is chosen asan important and tangible example of obtaining clear images through ageneral field variant distorting media. It should be evident that otherphysical instruments, subject to distorting mechanisms other than anatmosphere, can also produce sets of digital images and correspondingsets of field variant distortion estimates, thereby being capable ofutilizing the claimed methods of processing these sets of digital datainto a single clear image.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 pictorially summarizes some of the principal components of asystem according to a preferred embodiment of the present invention.

FIG. 2 is a schematic overview of the physical components of thepreferred embodiment.

FIG. 3 is a graphic representation of some of the mathematicalfoundations which are used to generate specific implementations of theprocessing methodology.

FIG. 4 is a detailed schematic of camera system embodiment using asimple wide field Hartmann camera embodiment.

FIG. 4A shows a typical array of sub-apertures placed at the prime focalplane of the FIG. 4 embodiment.

FIG. 4B is a schematic representation of spot images as they mightappear on the Hartmann detector in the FIG. 4 embodiment, with a furtherblow-up to a depiction of the centroids of those spots during typicalturbulent conditions.

FIG. 5 is a schematic functional drawing of a wide field wavefrontsensing system according to one embodiment of the invention.

FIG. 6 is a schematic representation of a field angle multiplexing beamsplitter.

FIG. 6A is a schematic view of the mirrors in the mirror sub-unit ofFIG. 6 as seen by the incoming light beam.

FIG. 6B is a front view of the two checkered mirrors of the FIG. 6embodiment.

FIG. 6C is a schematic drawing of the field stop sub-apertures which canbe used in a 2X wide field Hartmann embodiment.

FIG. 7 is a schematic view of a wide field Hartmann camera according toone embodiment of the invention, called the "full" wide field Hartmanncamera.

FIG. 8 is a schematic drawing of the optical train of one embodiment ofthe invention, along with a list of example specifications which areused throughout the disclosure in order to further explain detaileddesign issues.

FIG. 8A is a schematic drawing of the 16 round pupil images as theywould appear at the lenslet plane in front of the wavefront detector,with the lateral displacements of the beamsplitters overlaid.

FIG. 9 is a hardware and computer system block diagram of anillustrative embodiment.

FIG. 10 is an overview of the processing methodology used in anillustrative embodiment.

FIG. 11 is a representation of a generic weighted average operation.

FIG. 11A is a logical step overview of the role that the weightedaveraging operation can play in mitigating the traditional difficultiesinvolved with image estimation based on inversions.

FIG. 12 is a schematic diagram of the centers of field angles at whichthe wavefront is measured, overlaid on a grid of solution pixels.

FIG. 12A are examples of spectral sensitivity curves of the light beamswhich enter the primary camera (on the left) and the wavefront camera(on the right).

FIG. 13 is a schematic representation of the illustrative embodiment ofthe overdetermined basis functions used in the processing methodology.

FIG. 13A is a representation of how these basis functions are arrangedinside computer memory in an illustrative embodiment.

FIG. 14 is a diagram showing a processing methodology according to anillustrative embodiment of the invention.

FIG. 15A is an assistance schematic, clarifying some of the processsteps of FIG. 14.

FIG. 15B performs a similar clarification on another programming item.

FIG. 16 displays some of the general characteristics of "the borderproblem" referenced herein.

FIG. 17 is a processing step overview of illustrative methods whichdetermine image plane distortions based on Hartmann spot data.

FIG. 18 is a schematic which takes a closer look at the image dataderived from the full wide field Hartmann camera, showing two levels ofzooming in, culminating in a general look at nominal spot brightnessoverlap and crosstalk.

FIG. 18A looks at "in operation," i.e. turbulent conditions, spotcrosstalk and overlap.

FIG. 19 is a detailed description of the processing step of the globalwavefront sub-routine, which is a sub-routine in the main program of theillustrative embodiment.

FIG. 19A is a description of the processing steps involved in a localcentroding sub-routine.

FIG. 19B is a description of a sub-routine which estimates the pointspread function from a complex pupil function.

FIG. 20 is an intuitive graph which illustrates the relation between"obvious" applications of prior art methods (1A) to six novel methods ofthis invention (1,2,3,4,1B,&1C).

FIG. 21 is detailed description of the processing steps of the pointspread function method of creating "rev₋₋ array," an array used in themain program.

FIG. 21A is a description of the local wavefront sub-routine,corresponding to the point spread function method.

FIG. 22 is a graphical depiction of the relationship of novel processingmethods of this invention to the "obvious" application of prior artprocessing methods.

FIG. 23 is detailed description of the processing steps of the datapoint window function method of creating "rev₋₋ array," an array used inthe main program.

FIG. 23A is a description of the local wavefront sub-routine,corresponding to the data point window function method.

FIG. 24A is a schematic diagram assisting in the explanation of the datapoint window method of FIG. 23, depicting in particular the mirroring ofindices between psf array and the window array, in the lower right twodrawings.

FIG. 24B graphically provides further insight into the operation of theprocessing steps of FIG. 23, shown as they might affect the H matrix inequation 1 of the text.

FIG. 25 is a schematic with extended explanatory text in the figure,outlining the rationales of refining the method of determining rev₋₋array for use in the main program reversion process.

FIG. 26A is a schematic which highlights optical spatial filteringeffects to which the wide field wavefront camera described herein willgenerally be prone to. The lower figures detail the effective fourieroptics operation performed by the field sampling aperture.

FIG. 26B schematically outlines the framework in which the effects ofthe spatial filtering are mitigated by the processing methodology ofFIG. 27.

FIG. 27 is a preferred embodiment of the processing methodology which anmitigate the errors that would otherwise have been introduced by thespatial filtering of the sampling apertures of FIGS. 26A and 26B.

FIG. 28 is a general overview of other refinements to the dataprocessing.

FIGS. 29A-D show how geometric distortion manifests itself on a typicalgrid image, and how that distortion can be described graphically.

FIG. 30 is a simple example of how a 2 by 2 lenslet array might image agrid in typical atmospheric distortion.

FIG. 31 is a process flow chart for largely removing the effects ofscintillation in a distorted image.

FIGS. 32A-D show how one dimensional distortion can be measured.

FIGS. 33A-C comprise a graphic description of how two dimensionaldistortion can be measured using the "fat-to-thin" iterative lineintegral method.

FIGS. 34A-B show how a sparse field of brightness distributions mightappear, and how both photon flux and object contrast contribute to themeasurement of distortion.

FIGS. 35A-C are graphs which help explain one mode of finding thegeometrical distortion between a given distorted image and a referenceimage.

FIGS. 36A-B are graphics explaining the basics of how a single detectorsystem can manifest a changing distortion and focus on a given object.

FIG. 37 is the process flow of removing the higher order distortioneffects of a single detector system.

DETAILED DESCRIPTION

The disclosure begins with a general description of the aggregate ofparts of the system and then will discuss each part in detail. FIG. 1contains a schematic summary of the central aspects of the preferredembodiment. The physical entities are enumerated and the process flowwill be described in the text. In FIG. 1, an object 20 is being imagedby a telescope 6 with an atmosphere 3 distorting the light wavefronts 18emanating from the object, where the distortion on the wavefrontsappreciably vary depending on which point of the object the wavefrontoriginates. This variances according to field angle has been termed"field variance" and lies at the core of the novelty of certain aspectsof this invention relative to the prior art. The current inventionrecognizes that the incoming light wavefronts can be seen as arrivingfrom a grid of field angles 1--or more generally, from a functionalsteradian space. The grid spacing, or equivalently, the appreciablefrequencies of the field variance, is intimately connected to thestatistical characteristics of whatever atmospheric conditions prevail.The preferred embodiment of the invention therefore will explicitlymeasure and/or estimate the wavefront distortion as a function of fieldangle, using this estimate in subsequent data processing methodologieswhich create a clear image 5 of the object 20. The preferred embodimentdoes this by first utilizing the common methods used in adaptive opticssystems and speckle systems incorporating paraxial wavefront sensors.The primary beam of the telescope 6 is split into two beams by adichroic 28, where a narrow visual band beam falls on a primary detector8, while a generally wide band light beam enters a wavefront sensingsub-system, called a "wide field wavefront camera" in this disclosure inorder to summarize the overall difference from the prior art. The lightbeam which is split by the dichroic 28 and sent toward the wide fieldwavefront camera first encounters a field angle modulation unit 46.Several variants of this unit will be disclosed. The basic function ofthis unit is to effectively separate the identities of the various fieldangle directions, though the term "modulate" is the more appropriateterm since later image processing methods better view this function inthat mathematical context. A field angle modulated wavefront 48 emergeswhich then falls on a standard lenslet array 52, producing what iscommonly known as Hartmann spot images on a wavefront array 10. Theoutputs of the arrays 8 and 10 are digitized and sent to image storageareas labelled primary and wavefront in the figure. The wavefront datais first processed into an estimate of the field variant complexwavefront which exists at the pupil plane of the telescope. From there,it is further processed into a field variant focal plane estimate,stored either as a point spread function, itself a function of thewavefront field angle grid 1, or as its UV plane conjugate. Much of whathas taken place here utilizes basic techniques of prior art systems, butfurther involves the field angle modulation unit 46 and the use of the"FV" in front of the "pupil plane" and "focal plane," indicating novelprocessing methodologies which explicitly deal in field variantframeworks. A further novel aspect of this embodiment is then depictedinside the dashed box, where reference is made to the use of signal tonoise segregating basis functions. As will be seen in the more detaileddescriptions of the disclosure, a clear image 5 will be constructed byusing a generally overdetermined set of basis functions whichnevertheless provide for a common ground for individual "speckleexposures" to be compared and processed together in a weighted averagesense. It is interesting to note that these basis functions can reduceto the "trivial" case of the global fourier coefficients in imagingsystems which are "field invariant," whereas much more involvedengineering considerations lead to their specification in field variantsystems such as that of the invention. The dashed box in FIG. 1 thusdescribes how some set of E exposures, typically 1 to 50 for brighterobjects, are first transformed into these new basis functions, the fieldvariant focal plane distortions are likewise transformed into this sameset of basis functions (with certain extra considerations as thedetailed disclosure will show), a weighted average on the individualcomponents is derived which ultimately leads to a summed denominator anda summed vector product in the numerator (a weighted average operationgenerally has this form). After E exposures are cycled through, thenumerator is divided by the denominator--over all basis functions--andout comes the clear image 5 of the object 20. The importance of theseelements and processing steps drawn in FIG. 1 will become more clear asthe disclosure moves into more detailed levels of description.

FIG. 2 has an overview of the physical systems and sub-systems of thepreferred embodiment. It is perhaps most convenient to split thephysical components of the system into those elements which attachdirectly to the telescope and those which are part of an operatingstation. This disclosure concisely refers to this distinction as thecamera system 2 and the computer/data-collection system 4. The camerasystem 2 is attached to a large aperture telescope 6 with two electroniccameras inside, the primary camera 9 and the wide field wavefront camera11. These two cameras, in conjunction with camera controllers 12 andimage digitizing systems (framegrabbers) 14, are capable of collecting aseries of synchronized short exposure "digital images." Anywhere fromone to several thousand exposures will be gathered depending on a widerange of variables, where the larger numbers of digital images are oftencompressed and sent to mass storage 16, or otherwise processed "on thefly" if that is computationally economical. The exposure times (thelength of time the "shutter is open") for both cameras (which will besynchronized, incidentally) are typically on the order of one thousandthof a second up to as long as a tenth of a second or longer still,depending on many variables. It is well known in telescope design theorythat the light waves (18, FIG. 1) are collected by a large telescopeprimary mirror 22, reflected to a secondary mirror 24, and exit the maintelescope assembly as a panchromatic light beam converging toward anominal focal plane (exit beam not drawn). This emerging beam is thensplit by a dichroic 28 which passes a narrow spectral band beam into theprimary camera 9 and a broader band beam into the wide field wavefrontcamera 11 (further details will be described later, see FIG. 12A).During a given exposure, the primary camera 9 in conjunction with thecamera controller 12 and the frame grabbers 14 digitizes a wide field`distorted` image of the object in the narrow band of interest(typically 10-50 nonometers wide), while the wide field wavefront camera11 is explicitly designed, as will be described in the preferredembodiment of this sub-system, to gather complex phase information aboutthe effective pupil plane 34 across an array of field angles, not merelyalong the paraxial field angle. This wavefront beam is first measured bythe electronic array (10, FIG. 1) and converted to a digital image justas is the case for the primary camera previously described. As is wellknown in the art of wavefront sensing, the optics in the wavefrontcamera 11 are designed with the location of the effective pupil plane inmind (i.e. conjugated to the lenslet array 52, FIG. 1). Rounding out theitems in FIG. 2 are the cables which electronically attach the camerasystem 2 to the computer system 4, and then the computer workstation 15which essentially oversees and controls the operation of the entiresystem, as well as performing all or some of the data processingmethodologies described later on. The commercially available options forconfiguring items 12, 14, 15, and 16 are vast, and thus the blockdiagrams of these items are left in a general form. Typical componentswhich are implicitly contained in the computer workstation 15 areviewing monitor(s), disk drives, links to array processors orsupercomputers if needed, keyboards and other input devices, printers,etc. etc. Finally, the generic element of software 19 is depicted,encompassing both commercially available user programs, compilers, dataanalysis packages, AND, programs which embody aspects of this inventionwhich may either be on their own or incorporated into larger systemoperating software platforms (all as is well known in the art ofsoftware engineering).

As mentioned in the discussion surrounding FIG. 1, the physical systemin FIG. 2 will collect "E speckle exposures" in order to obtain oneclear image 5. E was mentioned as being, typically, anywhere from 1 to50 for brighter objects. An "exposure" in this case refers now to twodistinct synchronized digital images, one each deriving from the primarycamera 9 and the wide field wavefront camera 11. The pixel size or pixelextent of the two cameras need not be the same. Using currenttechnology, it would not be untypical that these exposures would begathered at a rate of one per second, thus the data gathered for oneclear final image 5 will take anywhere from a second to a few minutes.For applications where an object changes its intrinsic appearance at atime rate comparable to this rate or faster than this rate, specialfaster image capture systems would be needed in order to avoid timebased blurring of the final clear image 5. For a given exposure, thedigital image information gathered by the wide field wavefront camera 11is first processed by the computer workstation 15 into a field variantestimate of the complex phase of the light waves 18, FIG. 1 as they passthe pupil plane of the telescope 34; this is to say that the "fieldvariant pupil function" is being measured--precisely as conceived inprior art systems, only now it is explicitly a function of field angleas well. The field variant estimate of the complex pupil function isthen combined with a priori knowledge of the telescope optics, such asthe boundary of the pupil (its "support." including annularconsiderations), the magnification ratios applicable at the primarydetector array 8, and also possibly combined with knowledge of themicro-spatial sensitivity profile of pixels on the electronic array 8,and together a field variant estimate of the data point window functionis derived. The data point window function is the weighted field of viewthat each data point, or "pixel," has on the object space. Thesignificance of the data point window function and how it is derivedwill be explained in the processing section, but briefly, it is afunction which is in general much more useful in the correction of fieldvariant imaging problems (as opposed to using the point spread function)and it will play a central role in the data reduction steps of thesystem. A simple schematic of the data point window function is depictedin FIG. 3, outlining certain or its properties. FIG. 3 in many waysepitomizes a generalized digital imaging system, only that here theadditional window function modifier of the atmosphere greatly broadensand complicates the resulting window functions. The elements of thisdrawing begin with an index grid of data points 1000, where this is ingeneral a non-metric entity even though so many applications treat it ina metric fashion; the grid is of intuitional and practical convenience,corresponding to the nominal locations of pixels on a detector array.Each data point has a weighted spectral sensitivity function 1002 on ametricized object space 1004. The shapes and positions of these windowfunctions are determined by some set of modifiers or optical transferfunctions 1012, where the major ones listed here labelled are theatmosphere 1006, the overall instrument optics 1008, and the physicalpixel spectral/spatial sensitivity function 1010. Two particular windowfunctions corresponding to two random data points (pixels) are drawnwith projections onto object space which cross each other in deferenceto the inverting properties of typical imaging optics. The object planeitself is drawn as a square region corresponding to the nominalgeometrical projection of the primary detector 8 onto object space. Thisis purely as an intuitional aid as well. One subtlety alluded to by oneof the data points in FIG. 3 will be given more detailed considerationlater, and that is that data points near the edge of the detector havetheir window functions extend beyond the geometric projection of thedetector. Another item to note is that the window functions arethemselves smaller square regions than the entire object space. This isin reference to a technique for economizing on computational resources,whereby only a smaller region which contains some large amount--say99%--of the integrated volume of the function is utilized in processingstages. For the purposes of this general summary, it is enough tounderstand that the primary goal of the specific preferred embodiment ofthe invention is to measure these field variant window functions foreach exposure, and to use them at the multiplication stage of the dashedbox in FIG. 1 as well as directly in the calcuation of the denominatorin that dashed box. Highly refined procedures which will be outlined goeven a few steps beyond this.

In order to better understand what signal to noise segregating basisfunctions are and why they are used, we can briefly skip ahead to FIG.13 which has a simplified schematic depicted of the overdetermined basisfunctions used in the preferred embodiment, with one in particular blownup. This figure shows that the final image estimate is broken down intoa generally dependent (overdetermined) set of basis functions whichmaintain generally "knowable" signal to noise characteristics. Singularvalue decomposition techniques perform a more precise method ofsegregating the image into an independent set of basis functions, uniqueto each atmospherically distorted wavefront; whereas here, we merelywant to approximate this segregation in a simpler, computationallyefficient manner which is not unique to each exposure, hence the use ofoverdetermination. The purpose of this step is so that each individualprimary image exposure can separate out its unique signal to noisecomponents into a common framework so that it can add its own properlyweighted estimate to an averaged final image estimate. The basisfunctions chosen for the preferred embodiment are, effectively, theglobal fourier components AS CONTRIBUTED BY local data regions.[Computational efficiency will show that we don't need to use absoluteglobal fourier coefficients, however, hence only the 256 by 256 localfourier components]. In this manner, a given exposure with a given localregion which has a global "null" frequency will add no weight to a finalimage estimate, while another local region with a good signal to noiseratio on that frequency will add a proportionally high weight to a finalobject estimate. This separation and weighting of signal to noiseproperties is the essence of utilizing signal to noise segregating basisand has long been practiced in the far more intensive SVD matrixmethods.

FIG. 13 shows the primary detector plane 8 segmented into small squareregions of pixels which approximate a field invariant region, where onesuch region is projected outward. The basis functions are then definedin object space and wind up being the much wider region of UV planecoefficients (fourier coefficients) that would result as if the smallersquare segments were "zero padded" well beyond their actual extent allthe way to the 256 by 256 borders. This effective zero padding canextend as far as the extent of the entire object estimate if necessary,though a smaller region generally suffices for much the same reasonsthat the window functions 1002 FIG. 3 are smaller. This choice of basisfunctions effectively represents "the contribution of this smallquasi-isoplanatic segment" to the overall object estimate. Engineeringrefinement details can lead to a minimization of the errors inherent inthe quasi-isoplanatic approximation to any level required by overallsystem specifications: that is to say, the data sub-regions can get assmall as one pixel if the error criteria demand it and the computationaloverhead is acceptable. The choice of an overdetermined set of basisfunctions is a convenience and is generally dictated by the fact thatreasonable levels of turbulence lead to an ever increasingly difficulttask of using a "full independent" set of basis functions and stilladequately separating the signal to noise levels for generic atmosphericconditions. This just illustrates the point that normal engineeringconsiderations and refinement drive the system by system choices of whatexact basis functions to use, and that the exact choice of basisfunctions themselves is not a central aspect of the invention. Therecent article "Multidimensional digital image representation usinggeneralized Kaiser-Bessel window functions" written by Rober Lewitt,Journal of the Optical Society of America A, vol 7, No. 10 is just oneof many published articles which illustrate both the wide range ofchoices of inherent basis functions which could be utilized, AND, thedegree of analysis which can go into choosing a proper set.[Unfortunately, the window functions in the title are not at all thesame as the ones used in this paper. This will be discussed later, andhas no bearing on the point made here.] Though this article isostensibly written for a much different image processing problem, itnevertheless is representative of the extent and quality of refinementthat can be applied toward the use of basis functions in the service ofthis invention. In plain English, there is certainly room forrefinement.

Continuing with FIG. 13, each exposure is in turn broken down into thesedata segments, the segments are stepped through, the frequency domainrepresentation of the local window function--as determined by theprocessing of the wide field wavefront sensor data--is multiplied by thefrequency domain representation of the data in the segment, and theresult is accumulated into three image estimate data planes: the realnumerator, the imaginary numerator, and the magnitude denominator. [FIG.1 treats the complex numerator as one entity.] FIGS. 11 and 11Aconcisely describe the rationale behind these data planes, showing themto be complex components of a typical weighted averaging process. Themore detailed discussion surrounding FIGS. 11 and 11A will show that theadverse theoretical and practical effects of null components, which arequite well known in the art of image restoration, is entirely removed atthis point due to the cancellation of potential divisions by zero by anequal multiplication of that zero in the service of the weightedaverage. In brute force deconvolution terms, the acquired signalcomponent would require a division by the magnitude of the transformingcomponent (which could easily be a zero, see step 2 FIG. 11A). Whereasthe weighted average operation requires a weighting entity which isproportional to the square of the signal to noise ratio, which nicelyworks out to be approximately the square of the magnitude of thetransforming component as well in a typical white noise, uniformstatistics imaging systems. The remainder of FIG. 11A merely outlineshow to square the magnitude of the division in order to exactly cancelthe squared magnitude of the multiplication. It is also important tonote that this same result can happen with other basis functions where adivision by an attenuated component (basis function) can be preciselyoffset by the signal to noise metric of that very same in the service ofthe weighted average; it is not unique to fourier formulations of theproblem. In many ways, this whole weighted average process is a veryfortuitous relationship, in that the very randomness of the atmosphereand the concomitant roaming of nulls in the localized UV plane insuresthat all local UV plane points eventually get a healthy batch ofreasonable signal to noise levels which can readily translate into aweighted average. [Whereas static systems can have nulls stuck at thesame place no matter how many data realizations are made].

Finishing up with FIG. 13 and FIG. 1, then, all exposures are processedin this manner, an overall image estimate is represented as theaccumulated real and imaginary numerators and magnitude denominators,each segment is then stepped through where the division of the complexnumerators by the single denominator takes place to give a local UVplane estimate, then this estimate is either directly added to theaccumulating final image UV plane estimate, or alternatively, is fouriertransformed into the spatial domain and added into the accumulatingfinal object estimate in that domain. After all segments are cycledthrough, the result is a clear image 5. FIG. 1 with a Strehl ratioapproaching 1.0 and a signal to noise ratio determined by a wide rangeof variables including all the aforementioned design issues. The overallStrehl ratio performance is also a function of which design trade-offsare chosen and what the prevailing seeing conditions are.

Now that the general description of an illustrative system is completed,it will be more readily apparent how and why the sub-components functionas they do. These items will now be described in detail.

THE CAMERA SYSTEM 2, FIG. 1

FIG. 4 depicts a more detailed optical layout of the camera system 2.The primary panchromatic light wave emanating from the telescope optics26 first encounters an (interchangeable) dichroic, 28, whose normal istilted at a 45 degree angle relative to the paraxial ray. The dichroicis designed such that the close neighborhood of incidence angles about45 degrees will transmit light waves within a narrow band of wavelengths(typically 10 nanometers to 50 nanometers) toward the primary camera 9,while the complement of this narrow band (broader band light with theprimary narrow band generally excluded) is reflected along theperpendicular beam which is directed toward the wide field wavefrontcamera 11. It is optional to put secondary filtering elements, 36 & 38,in the path of either or both of the two new light beams, giving rise tothe properly spectrally-shaped light beams 30 and 32, if thecharacteristic spectral curves of the dichroic alone is deemedsub-optimal for any reasons. Requiring one 45 degree dichroic to filterboth beams to the desired spectral shape might be too much to ask, thusrequiring further spectral filtering of the individual beams. As anengineering point, the primary light beam 30 should be a narrow enoughpassband such that the image data resulting from the primary camera 9maintains reasonable signal levels at the higher spatial frequencies fora given level of ambient turbulence; while the spectral passbandcharacteristics of the wavefront light beam 32 going into the wavefrontcamera depends ultimately on the kind of wavefront sensing camera used,and in the case of a wide field Hartmann type used in the preferredembodiment of this disclosure, the spectral passband is tied to overallsystem performance specifications and residual error criteria in thewavefront reconstruction processing. But in general, the spectral shapeof the wavefront light beam 32 is a function of the correlation betweenthe phase distortion at these wavelengths to the wavelengths utilized bythe primary narrow band camera 9. This is considered ordinaryengineering practice even in prior art systems, i.e. the determinationof precise spectral characteristics of the lightwaves sent to a Hartmannsensor. It will also be further discussed in later sections.

The primary camera light beam 30 converges toward the primary electronicfocal plane array 8 as in any ordinary imaging system. Before reachingthe primary detector, the primary light beam 30 may optionally encountermagnification optics 58 which is often used to match plate scales toeconomical designs of the primary detector 8, such as pixel size and/orextent. A shuttering device 40 is placed in front of the primaryelectronic array 8. This shuttering device is synchronized (opens andcloses precisely coordinated with) the shuttering device on the widefield wavefront sensor 42. Both shuttering devices are optimized forvariable short exposure times, capable of operation over a very largeamount of shutter cycles (10's or 100's of thousands of cycles or more).A mechanical shutter remains as a viable choice for these two devices,but certainly anything else can be used, such as electro-optic shuttersor detector array electronic shuttering. The primary electronic array 8should also be grossly interchangeable if either a very wide range ofVIS/IR imaging is contemplated, or optimized quantum efficiency focalplane arrays are matched to narrower bands of the spectrum--e.g., a newfocal plane every 250 nm with its QE appropriately optimized. Theprimary electronic array 8 integrates the incoming lightwaves as anytypical electronic camera, converting light intensity patterns as itssurface into closely corresponding electrical signal patterns(capacitively coupled electrons typically). The description of the dataacquisition electronics and the camera controlling hardware will explainhow these electronic signals are converted into a stored digital image.

The preferred embodiment of the wide field wavefront camera is what canbe called a "basic wide field Hartmann camera" 44. It is enumerateddifferently here (not "11" as in FIG. 1) to denote a particular design,the Hartmann based design with the simplest field modulatingconfiguration, in the generic class of wide field wavefront cameras thatcould actually or conceivably measure the pupil function 34, FIG. 2 as afunction of field angle, be these "other types" previously invented orotherwise. In this "simplest configuration" design, the one which willsuffice for bright extended objects such as planets, the wavefront beam32--optionally spectrally filtered by filter 38 as mentioned--thenconverges toward a nominal prime focal plane 66 at which an array ofsub-apertures 46 is placed. Before reaching the nominal prime focalplane, however, it is again optional to place magnification optics 60 inthe path of beam in order to increase or decrease the effective platescale at the nominal prime focal plane 66 (as is common in mostinstruments which attach to generic large telescopes). These optics 60can be reflective, refractive, or some combination of both. Theygenerally allow for matching later optical scales to convenient andeconomical construction of other elements in the wavefront camera as awhole. FIG. 4A then depicts a head on view of a typical arrangement ofthe array of sub-apertures 46. The sub-apertures themselves 48, which inthis example are 3 arc second diameter circles with a square arrayspacing of 12 arc seconds, giving a total of a 4 by 4 square array ofthese sub-apertures. The spacings, the sub-aperture patterns, and eventhe opacity (which in FIG. 4A is binary in that it either completelytransmits or completely absorbs), are all variables which can beoptimized for particular telescopes, on particular nights (seeingconditions), imaging in particular spectral bands, and utilizingparticular data processing methodologies. It is recommended that thespacings chosen in FIG. 4A are not used uncritically, and the laterdescription on the "full field" embodiment of the invention will go intomuch more detail on determining preferred design specifications. Thearrangement depicted in FIG. 4A is generally adequate for seeingconditions which practitioners refer to as an r_(o) of 20 cm or better,and marginally adequate as seeing conditions worsen. A large part of thereason for creating advanced designs of the wide field Hartmann camerais to address much worse seeing conditions. The description of theprocessing methodology which reduces the Hartmann data will betterexplain that there is no absolute cut-off point on choosing thesevariables for given seeing conditions, and that it is merely a matter ofa increasing error and noise performance as the seeing conditions worsenbeyond the designed optimum. A "field angle modulated" light beam 48emerges from the sub-aperture array 46. The sub-aperture array itself issurrounded by a (generally multi-element) field lens 50 which forms animage of the telescope pupil plane 34 onto a plane at which a lensletarray 52 is placed. The field lens 50 is drawn as a double elementrefractive type only symbolically. It is well known in opticalengineering that both reflective, refractive, and multi-elementsolutions to this field lens can be utilized. [Conceivably, this couldbe a single element of etched reflective circles deposited on an offaxis paraboloid surface having a background absorbing material; etc.Also, there appears to at least be the option of using no field lenswhatsoever if for some reason that might be preferable. It could save onreflective light losses, but the consequent slight blurring of the pupilimage may ruin wavefront reconstruction accuracies. This choice isessentially using the sub-apertures 48 as pinhole cameras. It is anoption, anyhow.] The lenslet array 52 sits at what is commonly referredto as the pupil image plane, i.e. the focal plane corresponding to theeffective pupil plane 34 as imaged by the field lens 50. The lensletarray 52 is typical of the kind which are readily availablecommercially. It may be noted that the arrangement of this lenslet arrayUSUALLY SHOULD BE SQUARE rather than hexagonal or anything else. This isbecause each lens in the lenslet array is imaging a 4 by 4 array ofsub-apertures 46, and the square packing of the lenslets would nicelyfit the square pattern of "spot images" resulting on the wavefrontdetector array 10 (see FIG. 4B). Since this is a wide field device whichis imaging, in this example, a 4 by 4 array of sub-apertures 46, thedesigned position of the lenslet array relative to the field lens 50(i.e. the field lens focal length), the aperture size of each and everylens in the lenslet array, and the f number of each and every lens inthe lenslet array must be chosen such that each lens in the lensletforms a non-overlapped image of the 4 by 4 sub-aperture array 46[overlapping images of the 4 by 4 sub-aperture array 46 are possible,but are not part of the preferred embodiment] Each lens in the lensletarray 52, where typically there are hundreds of lenses if not thousands,forms its own unique image of the sub-aperture array 46 onto thewavefront detector array 10. It can be readily seen that this willamount to thousands of spots being imaged onto the array 10. Eachindividual spot is, by design, "nominally allocated" an area on thedetector array 10 anywhere from a 2 by 2 up to a 10 by 10 array of"pixels," enough so that in-operator spot wander and normal spot spreadare adequately accounted for. Advanced engineering, where an intricaterelationship between the physical design and the processing methodologyis utilized in order to attain low residual wavefront error performance,will show that this idea or "pixel allocation" to each spot is somewhatoversimplified but is nevertheless quite helpful in the early design offull systems. FIG. 4B depicts a simplified view of the multitude of spotimages 54 as overlaid on the border outline of the wavefront detectorarray 10. FIG. 4B then zooms in on a small neighborhood of spots wherethe in-operation wandering centers of gravity, relative to the nominaldistortionless centers of gravity (i.e., the "null positions"), aredepicted as x's and points, respectively, on a pixel grid. [This is allwell known is the art, with the twist that these spots are actuallygrouped into individual images of the 4 by 4 sub-aperture array 46.] Thelarge circular pattern of spots depicted in the upper left of FIG. 4B isa simplified idealization. In actuality, these spots are noise pronephoton counts in an array of pixels with inevitable cross talk betweenspots. Assuming a field lens 50 is used, the spots would not have suchan abrupt circular edge but would instead "dim" together as square unitsfor the 4 by 4 sub-images near the edge. There would be an inner circleof no spot images also, if the primary/secondary mirror configuration ofFIG. 2 is used. Finally, the 4 by 4 sub-images created by the individuallenses in the lenslet array 52 have an exaggeratedly sparse packingdensity in order to graphically demonstrate their presence. Actualdesigns would pack these sub-groups much closer together to the pointwhere the individual 4 by 4 images become virtually indistinguishableand it becomes one large array of spots. It is an issue of refineddesign (which will be further discussed in the processing descriptions)to determine at what separation the spots must be to achieve a crosstalk level which the overall system specification can tolerate.

Referring back to FIG. 4, the shuttering device 42, previouslymentioned, could conceivably be placed anywhere in the Hartmann opticaltrain and in this example is placed in front of the lenslet array 52.Its chief function is to ensure that the collected optical wavefrontdata corresponds exactly in time to the primary camera data. Once again,the wavefront detector array 10 converts the resultant intensity patternof spots into a closely corresponding pattern of electronic signalswhich the data acquisition hardware will convert into a digital image.

Now that a simple but slightly limited form of a wide field wavefrontsensing system has been described in the context of FIGS. 4, 4A, and 4B,it will easier to explain the more fundamental aspects of a wide fieldwavefront sensing system. FIG. 5 discloses two primary elements of thewavefront sensing system of this embodiment of the invention whichdistinguish it from its paraxial (or mono-axial) prior art parents. InFIG. 5, it can be seen that a field angle modulation unit 47--where inFIG. 4 it wound up being the field lens 50/sub-aperture array 46combined unit--is placed (centered about) the nominal prime focal plane66 of the overall telescopic system. This unit 47 (or system, as in thecase of the preferred embodiment of this invention), explicitlymodulates the optical field existing at or near the prime focal plane 66in a manner that will allow later "traditional" wavefront sensingmethodologies to distinguish between field angles or the incomingoptical light waves. Moreover, the field angle modulation unit 47explicitly spatially filters, in the optical domain, the incomingwavefront light beam 32. This explicit spatial filtering is bestrecognized in classic fourier optical terms, being the transform planebetween the conjugated planes of the pupil plane 34 and the pupil imageplane at which the lenslet array 52 is placed in FIG. 4. The secondcentral element to the illustrated embodiment of this invention is theaddition of a processing domain stage 41 which explicitly recognizes thespatial filtering which is impressed on the optical wavefront by theunit 47, and applies a pseudo inversion of the filtering in order tobetter estimate the actual wavefront function existing at the pupilplane 34. It is possible to exclude this processing step 41 at theexpense of higher residual wavefront error. The second on processingmethodology will have a great deal more to say about box 41. Ancillaryto these two central elements 47 and 41, reference is made to "modified"traditional wavefront sensing systems 43, since the fact that the fieldangles are modulated and spatially filtered can lead to specificengineering design changes made to otherwise "traditional" prior artsystems. In the case of FIG. 4 and the simple configuration, thismodification was made to the traditional Hartmann sensing prior art,where it was necessary to explicitly design around the fact that it wasactually 4 by 4 sub-arrays of spots which were being imaged. This"modification" becomes more pronounced in the preferred embodiment ofthe wide field wavefront camera presented in FIG. 7. FIG. 5 therefore isa basic blueprint from which a wide variety of engineeringconfigurations can derive. Principal variables include: field anglemodulation methods, spatial filtering methods, beam splitting methods inthe preferable systems which attempt to create "lossless" measurements(see FIG. 7), the actual type of prior art wavefront sensing methodsused, and finally the type of processing inversions which can beapplied. The preferred embodiment of this invention has made quitespecific choices regarding these many variables, choices which attemptto balance a wide range of current economic considerations withperformance considerations. As with all engineering realizations of agiven type of invention, it should be expected that any given designwill have driven by its own unique set of performance and economicfactors. We therefore turn once again to the basic Hartmann sensorphilosophy of wavefront sensing systems for yet improved types of widefield wavefront systems, recognizing that there are quite viablealternatives that could derive from the embodiment depicted in FIG. 5.

A fundamental variant to the simple wide field Hartmann design of FIG. 4is what can be referred to as the beam splitting wide field Hartmanncamera. It is entirely probable that of the two specific beam splittingdesigns about to be described, the second or "full field" design will bethe most commonly utilized in physical systems: even though it is themost complicated and expensive, its performance gain will usually beworth it. Therefore, of these, the full field design is presentlypreferred, while the "simple" design presented in FIG. 4 and the "2X"design about to be presented will be viable alternatives to the fullfield design. (A third implementation is described at the end of thisspecification.) [The 2X design in particular may some day get very closeto the full field ideal, especially when yet further refinements aremade to "descrambling" the information gathered by the wavefrontdetector 10. More on this in the section describing the processingmethodology.] The simple configuration listed in FIG. 4 was describedboth as a modesty functional option, and, as an easy introduction to theprinciples involved. The primary considerations for considering the beamsplitting designs are: a) sensitivity, b) spatial filtering, c)guaranteed field coverage of bright (beacon) objects, and d) field anglesample spacing. The essential idea of the beam splitting variant designsis to multiplex the field angles out along separate beams so that eitherthe field angle sampling aperture diameter can widen, or the field anglesample spacing can decrease, or some traded off combination of both. Anaccompanying principle is to gather as many "signal photons" aspossible. The widening of the sub-aperture diameter effectivelyincreases the signal photons available, reduces the spatial filtering ofthe re-imaged pupil function at the lenslet plane 52, and helpsalleviate potential difficulties involved with bright stars (beacons)"falling between the cracks," i.e. falling outside the sampling windowof a given field angle. The ability to decrease the field angle samplespacing translates into the ability of the device and of the overallsystem to function at increasingly worse turbulence conditions (lowerr_(o) 's). This latter ability is arguably the largest value of thesebeam splitting variant designs.

The new element of a "2X" device, in relation to the simpleconfiguration of FIG. 4, is depicted in FIG. 6 where a new opticaldevice referred to as a field angle multiplexing beam splitter 56replaces the combined unit of the field lens 50 and the sub-aperturearray 46 depicted in FIG. 4. In effect, this new device 56 is thegeneric field angle modulation unit 47, FIG. 5. The positioning of thisnew device 56 is referenced by the relative position of the nominalprime focal plane 66, whereas detailed design should reference thenominal prime focal planes drawn in FIG. 6A, discussed later. The 2X ofits title refers to the fact that it can halve the sampled field anglespacing for a given sub-aperture diameter, relative to the simplestconfiguration of FIG. 4. The beam splitting device 56 in FIG. 6 has, asa central element, a carefully arranged series of 45 degree angled fullmirrors 74a-d and checkered mirrors 64 & 70 which, rigidly held togetherby a mechanical holding framework and collectively called the mirrorsub-unit 81, produce four outgoing beams each with a virtual 4 by 4sub-aperture array variously offset from each other, i.e. spatiallymultiplexed, giving a total of 64 sampled field angles. This arrangementas presented here allows for the same sized sub-aperture diameter to beused as in FIG. 4, i.e. a 3 arc second diameter circle, but now placedon 6 arc second spacings. This keeps the effective sensitivity and theeffective spatial filtering effects of each field angle sampleessentially the same, except now there are a few more reflectivesurfaces which will contribute to light loss. If need be, the aperturescould grow to, say, 5 arc seconds in diameter and the field spacing to10 arc seconds. Judging from typical anisoplanatism measurements andcalculations, a 6 arc second field angle spacing could perform areasonable field angle sampling down to r₀ =10 cm or perhaps lower,while the 10 arc second spacing would be adequate for proportionatelylow turbulence levels. If the device of FIG. 6 is precisely substitutedfor the field lens 50/sub-aperture array 46 of FIG. 4, then theresulting four light beams 75a-d emerging from the device likewise fallonto the lenslet array 52, with the remainder of the system being aspreviously described. Though FIGS. 6A-C will describe the mirrorconfiguration in more detail, FIG. 6 points out that there is a single"half" 51 of the field lens functional system for the incoming wavefrontlight beam 32, while there are four separate lens systems 53a-d for theoutgoing four light beam 75a-d. Again, the drawing uses single lenses assymbolic of the field lens function, while common optical principlesshould be applied to the detailed design of the overall field function,including chromatic aberration control as is common; furthermore,reflective elements are certainly feasible solutions as well. [In FIG.6, only three of the four outgoing lenses 53a-d are drawn forsimplicity, there other one being hidden] A final note on the fieldangle multiplexing beam splitter 56 is depicted in FIG. 6 by the plug 79and the encasement around the entire device. This represents that theentire device is sealed and filled with index matching fluid, matched tothe (special mean) index of refraction of the glass flats that are usedas bases for the checkered mirrors 64 and 70 (again, FIG. 6A-C will gointo more detail on these). The inner lens pieces of the field lenssystems will also have, in general, more closely related indices to thefluid than they would to air. The reason for the index matching fluid israther obvious: reduce unnecessary photon loss. The function of the beamsplitting mirrors is purely reflective, and thus reflective losses mustbe tolerated. Refractive losses on the other hand are completely a sideeffect of the need for a refractive substratum, and the index matchingfluid essentially removes this side effect. The beauty of the use of theindex matching fluid is that it all but offsets the overall photon lossthat the beam splitter 56 would have otherwise introduced relative tothe simple configuration. Removing the field lens air-glass glass-airinner surface reflective loss of the simple configuration is a partialcompensation for the 2 to 4 reflective surface losses experienced byeach of the four beams (one beam has two reflective surfaces, two havethree reflective surfaces, and one has four reflective surfaces).

Obviously there is much more to say about the mirror sub-unit 81 and itsfunction. FIG. 6A therefore has a top view of the detailed arrangementof the preferred embodiment of the mirrored surfaces, i.e., it is seenfrom the perspective of the incoming wavefront beam 32, whose centralray is depicted as the x on the first redirection mirror 62. There arecertainly alternative ways to arrange full mirrors and checkered mirrorsin order to split the incoming beam into four new beams, including theuse of four detector arrays for each of the four beams, but thepreferred embodiment attempts to make efficient use of a single detectorarray 10. For some applications on very large telescopes, however, theremay be no choice but to use multiple detectors merely because theremight not be any detector available with enough pixels to image all theresulting spots. The surfaces of both the redirecting full mirrors 62and 74a-d are typically aluminum (if using visible wavelength light)deposited or etched onto a thin precision glass flat using generallyavailable commercial processes. The optionalmagnification/demagnification lens 60, FIG. 4, can potentially create animage size of the 8 by 8 array of field angles which is optimized forease of manufacturing the glass pieces, depositing patterned reflectivesurfaces upon them, cutting them into the proper shape, and finallyplacing them in precision relative positions via an appropriatemechanical housing which stays clear of the light paths. For thepurposes of exposition, we can image that manufacturing considerationsof this sub-unit 81, FIG. 6, leads us to specify a virtual image size of1 square centimeter for the 8 by 8 array of field angles (i.e. 8×6" perangle giving 48" square steradian over 10 mm by 10 mm); the telescopeoptics (f number) and the size of the wide field wavefront electronicsensor 10, among other things, play a large role in best determining thevirtual image size at this stage as well. Referring again to FIG. 6A,the incoming light beam 32 is symbolized by an x in the center of thefirst redirection mirror 62, directed from above the plane of thedrawing toward the plane of the drawing, then reflecting off the firstfully mirrored surface 62 which is a front aluminized flat with itsnormal directed at 45 degrees to the plane of the drawing. Its mirroredsurface is on the side of the flat which is closest to the viewer, i.e.above the plane of the drawing (drawn as clear). The light wave 32 isdirected by the mirror 62 downward in the plane of the drawing, depictedby an arrow emanating from the x toward the checkered mirror 64. It thenencounters this first checker patterned mirrored surface 64, depicted infront view in FIG. 6B, left drawing. This surface is an 8 by 8 checkerpattern with one dimension being the square root of two longer than theother, giving a projected square pattern when seen from a 45 degreeangle. The darker regions are a thin layer of aluminum (or other chosenreflective material) which will reflect most of the incident light whilethe white regions are clear glass, allowing full transmission of thelight. Referring back to FIG. 6A, this mirror unit is placed at a 45degree angle immediately preceding the nominal focal plane 66, depictedin that figure as two separate nominal focal planes 66a&b due to theselective splitting of the checkered mirror 64. It can be readily seenthat the alternating reflecting regions and transmissive regions will beincreasingly displaced from the nominal focal plane as they extendoutward, i.e., the unit 64 is by necessity skewed and increasingly outof focus. In general, this displacement will take place within the usualand somewhat arbitrary "depth of field" about the nominal focal plane66. More precisely, however, the overall system effects of thisdisplacement is a rounding of the sampling window on the distortion of agiven field angle. As will be better explained in the description of theprocessing methodology, the adverse effects of this asymmetric blurringis virtually negligible. This can become more obvious if one imaginesthat we replace the simple checker pattern of FIG. 6B with aluminum"densities" which are sinusoidal in nature, giving rise to a spatiallyvariable beam splitting rather than the "hard edged" one employed in thepreferred embodiment. This option will indeed be explored in the courseof normal engineering refinement and is at the root of why the device 47in FIG. 5 is termed a "modulation" device as opposed to, say, a"splitting" device; modulating the light waves at the nominal focalplane is just as acceptable as "separating them"; the advantage of theapproach of the preferred embodiment is that many techniques can now beborrowed from prior art Hartmann systems. In this case, there is no needwhatsoever for "edges" and the displacement about the nominal focal 66is of less significance still (see FIG. 26B also). Besides, it is highlyrecommended--though not absolutely necessary--to place two field stops73a&b at the two respective nominal focal planes 66a&b; (in very lowturbulence situations it might be possible to remove these, seeprocessing methodology section). These field stops 73a&b, where one isdepicted in FIG. 6C while the other is just this same design onlyflipped about the horizontal axis, create a well known spatial lightwavefiltering which will later be used by the processing methodology anddata reduction. It also helps keep spot crosstalk under control. Twolight beams 68a&b (drawn as two arrows crossing through the nominalfocal planes 66a&b) now emerge from the first checkered mirror 64, oneimpinging on the left half of the double checkered mirror 70, and theother impinging on the right half (left and right referenced by FIG.6B). FIG. 6B also has a front view drawing of this double checkeredmirror 70. It can be seen in FIG. 6A that the two beams 68a&b arefurther split apart into 4 light beams enumerated collectively as 72a-d.A detailed look at the patterns will exhibit that they are offset fromeach other such that, for each contiguous 2 by 2 array of field angleswhich entered the system, it is now broken out so that each of the fourincoming field angle travels into only one of the four outgoing beams72a-d. Thus, the field angles have been spatially separated and thenmultiplexed. Overall, it can be seen that the new beams 72a-d consist ofa four by four array of field angles with twice the field angle spacingsas would be the case in the simple configuration of FIG. 4 (we couldhave widened the aperture diameters, 73a&b, as traded off by field anglespacing). It can be seen that the checkered mirror 64 is butted upagainst the center of the "double checkered" mirror 70 at a normal toits surface. The new beams 72a-d each encounter a 45 degree angled fullymirrored surfaces, called the 4 outgoing redirection mirrors 74a-d, withtheir mirrored surfaces on the other side of the viewer, thus directingthe new beams down into the plane of the paper, back into the originaldirection of the incoming beam 32. Finally, as an intuitive aid in howand why such a device as this functions properly, there are large dashedcircles 76a-d which depict the image of the telescope pupil function asprojected by these redirecting mirrors 74a-d onto the plane of thelenslet array 52 in FIG. 4 (i.e., as if the plane of the lenslet arrayis some distance into the plane of the paper). The lateral distance ofthe mirrors 74a&d from the central assembly is such that these circlesalmost touch (to within engineering tolerances); the circles on theupper left/lower right axis in FIG. 6A are further offset due to thefirst beam splitting. Good engineering practice would minimize thisoffset by keeping the plate scale at the nominal focal plane as small(least magnified) as is reasonable from a tolerance andmanufacturability standpoint. It can be seen that a square electronicfocal plane array 10 used behind the lenslet array 52 should be tiltedat a 45 degree angle relative to the axis of the incoming 8 by 8 arrayof field angles in the light beam 32, in order to efficiently containthe multitude of spot images which will roughly correspond to thepattern of the large circles 76a-d. A few final engineering notes are inorder. First, all of the mirrors 62, 64, 70, and 74(A-D) are held inplace by a precision optical mounting structure (part of the sub-unit81, FIG. 6A) which would need to steer clear of the light beams andtheir widening breadth as they get farther from the nominal focal plane66. No completely reliable rule of thumb exists which specifies howprecise, in terms of optical wavelength tolerance, these surfaces needto be, nor their relative alignments and relative tilts. It is fair tosay, however, that the flats upon which the mirror surfaces aredeposited will need to be of high precision optical grades which willnot add their own distortions beyond the overall system specificationfor the wavefront sensor. The tilt tolerances appear to be lesscritical, in general, but this statement is rather arbitrary itself. Thepoint is that this is an area for precision optical engineering,informed by the informational needs of the processing methodology andthe overall system error budget. One true subtlety which may or may notbe necessary to address relative to overall system design would be toplace a trapezoidal "fanning out" of the dimensions of the checkeredmirrors on 64 and 70 in order to precisely track the slight divergenceof the incoming light beam (according to the F number of the incomingbeam 32); for higher f number systems this may not be such a bad idea;in that case, the checker pattern would appear to be a trapezoid on itsside in FIG. 6B. As another engineering subtlety, the spacing betweenthe checkered surface on unit 64 and the left half checkered surface onunit 70 needs to be far enough so that the beam 72B does not clip theredirecting mirror 62. This is a particular instance of the general needto keep the mechanical surfaces out of the way of all the light beams.Finally, the slight bulge drawn in FIG. 6A on the checkered mirrors 64and 70 are an exaggerated thickness of the underlying metallicreflecting material. In fact, this is the recommended side, though oncethe consideration of index matching fluid is introduced, the choice ofwhich side of the glass flats should have the reflective materialbecomes a design option.

The next and most comprehensive variant on the nominal design of thewide field Hartmann camera depicted in FIG. 4 is what will be referredto as the "4X" wide field Hartmann camera, or better, as the fullcoverage wide field Hartmann camera, depicted in FIG. 7. The guidingprinciple behind this device is that no photons are lost other than bysimple reflection on refraction. The two previous configurations areboth forced to truncate (occlude) the field of view along the centralaxis of each field angle in order to maintain proper segregation of thespots on the electronic detector array 10. The full coverageconfiguration in FIG. 7 will be shown to multiplex the incoming beam 32into a 4 by 4 array of variously offset new beams--16 new beams in totalfalling on the lenslet array 52. As will be seen, this is a sufficientsegregation to allow for the sub-aperture breadths of each field anglesampling window to grow to the size of the field angle spacing itself,thereby allowing for a full coverage of the field of view (no brightstars will fall between the cracks; they may fall "on a crack," i.e. asomewhat sharp transition zone, but not between the cracks). Moreimportantly, it will allow for the maximum sample spacing density of thefield angles (relative to spot overlap and crosstalk limitingspecifications).

Referring now to FIG. 7, the full coverage configuration is identical tothe 2X configuration described previously, up through the 2X beammultiplexer 56, except that the sub-aperture field stops 73a&b FIG. Cwill not be placed at the nominal prime focal planes 66a&b. The Hartmannbeam 32 passes through the optional magnification stage 60(demagnification being an implied possibility as well). It enters andexits the beam multiplexer 56 just as described above. As mentioned, theone definite change to the interior design of the beam multiplexer 56 isthe removal of the field stops 73A&B. The full coverage device will alsorequire some adjustment to the lateral distances between the redirectionmirrors 74a-d, FIG. 6A, and the central assembly in that same figure.This readjustment in lateral distance has to do with the fact that eachof these new beams will again be split into 4 new beams each, giving atotal of a 4 by 4 array of new beams 88a-p, and so the nominal axiscenters of the outgoing beams 75a-d need to account for this. Gettingback to the first stage, the four light beams 75a-d, laterally displacednow, travel from the beam multiplexer 56 to the first pupil image plane90 where the lenslet array was placed in both the simple configurationand the 2X design. For the full design of FIG. 7, however, this firstpupil image plane 90 has four new field lenses 98a-d aligned with thefour beams 75a-d. These new field lenses can again be refractive orreflective as a designed may choose, where here they are depicted asreference. These lenses are drawn crudely in isometric view in order tosimplify their depiction and should be seen as centered about the sameplane 90. These lenses are also roughly drawn as a symmetric 1 to 1finite conjugate lens system relative to the first prime focal plane 66and the second prime focal plane 93, but if there is determined to be aneed for a remagnification at this stage due to engineeringconsiderations, it certainly could be incorporated in this lens system92a-d. Reiterating earlier comments, the design issue of magnifyingand/or demagnifying the various optical stages is driven the primarytelescope f number, the preferred sizes of the beam splitting mirror andtheir mechanical mounts (cost of manufacture), and the pixelsize/overall size of the focal plane array 10. These, and other normalengineering considerations, are all significant factors in determiningwhere and when to magnify or demagnify the succession of image planesand pupil planes. It should be noted that in consideration of minimizingreflective and scattering photon losses, these magnification stagesshould be kept to an absolute minimum. The reimaging lens system 92a-dreimages the nominal focal plane 66 onto a new nominal focal plane 98.About this new nominal focal plane 98 is placed an array of fourbeamsplitting multiplexers, 100a-d, of the same design type as 56, hereenumerated as 96 being the entire four unit array. The mechanicalpositioning and centering of the array 96 corresponds to the position ofthe outgoing beams 75a-d and the corresponding positions of their imagesat plane 98. There are only three key differences between the design ofone of the second stage beamsplitters 100 and the first stage beamsplitter: 1) (optional actually, if one to one mag not chosen by fieldlenses 98a-d) overall scale; 2) 4 by 4 checkered mirrors rather than 8by 8, giving rise to a doubling of the dimensions of the individual"checkers" on 64 and 70; and 3) lateral spacing of the four redirectingmirrors are now optimized for packing spots onto the array 10.Significantly, all four beamsplitting multiplexers 100a-d are identicalto each other in design and manufacture. Also, none of them containsfield stop arrays, 73(A&B). The only true subtlety involved in thenominal alignment of the array 96 with the incoming beams 74a-d is NOTto reference the absolute central axes of the incoming beams, butinstead to reference the center of the unique offset grids contained ineach beam. [Referring to FIG. 6B, the outgoing images are similar to themirror patterns on the double checkered mirror 70; the beamsplitters100a-d should each reference it own uniquely offset grid pattern] Themultiplexers 100a-d also contain index matching fluid, surrounded by afield lensing system which reimages the pupil image plane 90 onto thelenslet plane 52. There is yet another optional mag/demag lensing system(not depicted) possible between the beamsplitters 100a-d and the lenslet52, just in case the output from the array of beam multiplexers 96 hasan optical incompatible with the physical pixel sizes and extent of thefocal plane array 10. It is also conceivable that clever optical designcould turn the back element of the field lenses of the beam splittingmultiplexers 100(A-D) into a veritable mag/demag stage. This iscertainly an issue for normal detailed engineering refinement. Theremainder of the full coverage wide field Hartmann design is identicalto the nominal design in FIG. 4. It is up to the processing methodologyto sort out whose spots belong to whom, i.e. which beams are which andto what nominal field angle do they belong. As a final note to the fullcoverage design which has already been alluded to earlier, it may benecessary to consider putting an array of focal plane arrays, 4 focalplane arrays in total, one each at the output stages of each of thebeamsplitting multiplexers 100a-d. The potential need for thisconfiguration derives from the possibility that there are no acceptablylarge sized focal plane arrays available which can adequately containall of the spot information. The possibility even exists of placing anarray of 16 focal plane arrays behind each individual beam of theaggregate 16 beams 88a-p. These issues are exclusively in the realm ofdetailed system design and current technology, but the distinctpossibility to create these "arrays of focal plane arrays" exists. As aquick specific example of how a system might come to require such an"array of arrays," we can posit a 10 meter near full aperture telescopedesigned to work down to seeing conditions commonly known as r₀ =7 cm.System design parameters might likewise dictate the "allocation" of a 5by 5 pixel array per spot. Finally, there is an average 20% overallengineering tolerance beam separation between the pupil images of theoutgoing beams 88a-p as projected onto the lenslet plane 52; this Fig.accounts for several factors. Sraightforward engineering calculationsshow that there needs to be a linear diameter of 160 subaperture lensesacross the pupil image (10 m/7 cm), where each lens images a 2 by 2array of field angles giving a linear breadth of two spots. The totallinear pixel breadth per beam, 100(A-P), is therefore160×2×10(pixels/spot) or 3200. Multiply this by 4 linear breadth ofbeams across our hypothetical single detector, and this gives 12800pixels, then add 20% for engineering tolerances and this gives a 15360by 15360 pixel detector needed as the focal plane array 10. This is welloutside the reach of current state of the art focal plane arrays.However, a 4 by 4 array of 4K by 4K focal plane arrays is ostensiblywithin the reach of current state of the art, if not pocketbook. [Itshould be noted that this is an extremely rough engineering example heredescribed only in order to fully illustrate the possible need formultiple detectors. Even the orders of magnitude involved could besignificantly altered after a thorough system design has been completed:i.e. less engineering tolerance could be allowed through tighter design,or better processing methodology could further reduce the pixels perspot through clever information extraction procedures. Finally, thesystem could be designed toward a nominal r₀ =10 cm or more, and make dowith a gradually decreasing resolution performance (from neardiffraction limited) as the seeing worsens beyond that point. The ONLYpoint of this example is to see that the idea of using more than onedetector is not so hypothetical.]

For the remainder of this disclosure, the full wide field Hartmanndesign depicted in FIG. 7 will be assumed and utilized as the preferredembodiment, while the simple and 2X configurations will be mentioned assignificant alternatives with the same principles involved. This beingthe case, we turn now to a relatively detailed design of a full widefield Hartmann camera in order to illustrate the application of thedesign principles outlined above. FIG. 8 has a schematic breakout of thecritical specifications for the example telescopic system which will beused throughout the rest of this disclosure. This example design shouldnot at all be considered as an optimized design for a 2.5 metertelescope, it only uses simple numbers and fractional relationshipswhich can best illustrate the principles of design. Those practiced inthe art of telescope system design will understand that gross systemspecifications such as these are only the starting point for detailedengineering refinement. In fact, detailed investigation of thesespecifications will find that very little engineering tolerance has beenincluded. This was done to further understand the nominal relationshipsof the basic designs parameters. Also included in the example systemdesign is a few significant specific variations on the "symbolic"outline of FIG. 4 with FIG. 7's full design substituted. Thesevariations illustrate the points made earlier on how specific designscenarios (such as a 2.5 m telescope w/ 20 arc second square field ofview--FOV) can give rise to specific optical engineering designrealizations. This point is most evident in the realization of thedesign of the field lenses 51 and 92. In any event, the specificationsof FIG. 8 contain a sufficient set of information to illustrate theimportant aspects thereof. I turn now to discussing these specificationsin detail.

FIG. 8 is best described by starting with the specification of theprimary wavelength λ=500 mm. Those practiced in the art will recognizethat this is in actuality a specification of the central wavelength ofthe narrow spectral band light beam which will fall upon the primarydetector 8. Details of the spectral characteristics are not pertinent tothis discussion and are contained in several other places in thedisclosure. It should be understood, however, that typical physicalrealizations of imaging systems according to the invention will have aseries of selectable bands at which the system will operate, includingthe choices of "wide and narrow" bands centered about the samewavelengths. The changing of these bands may likewise change othercomponents of the system, even the primary detector 8 (as was previouslypointed out). The aperture shape and size of the example telescope willbe a circular aperture with an outside diameter of 2.5 meters. For thisdiscussion, it is immaterial whether this is annular as in FIG. 1 or afull aperture telescope. The inherent focal length of the primarytelescope is chosen to be 25 meters, giving an F/10 primary opticalsystem. The primary beam is split by the dichroic system 28, 36 and 38as described in the preceding disclosure. Looking first at the primarybeam, this example design does in fact choose to use the magnificationoptics 58, FIG. 4. These optics effectively change the incoming F/10beam into an F/40 beam which will fall upon a 1024 by 1024 square CCDarray 8. These optics, 58, can be reflective or refractive of course,and they can certainly be optimized for the current band being imaged,interchanged as the primary spectral bands are interchanged. Theindividual pixel size on the primary CCD 8 is chosen to be 10 microns.With a plate scale of 1"=0.5 mm, this gives an effective pixel spacingof 0.02" at our 500 nm central wavelength, a precise match to thediffraction limit Nyquist spacing. The detector as a whole coversapproximately a 20 arc-second by 20 arc-second square field of view;Since the wide field Hartmann camera about to be discussed has targetedseeing conditions of r₀ =7 cm, this field of view can be seen to be muchlarger than the isoplanatic patch size at this level of turbulence.[There is no question whatsoever that wider fields of view could bechosen or larger telescopes. This example was chosen because many of thespecifications fall out into nicely understandable forms; and theyconform to "obviously commercially available" system design components.]

Turning now to the full wide field Hartmann camera in FIG. 8, it hasbeen chosen that no magnification optics are used and the systemdirectly uses the F/10 of the primary telescope. This gives a platescale at the first primary focal plane 66 of 1"=0.125 mm. Working nowfrom the overall specifications listed below, we find that an 8 by 8array of field angles will be sampled with a 3.6" square spacing betweenthe field angles. This gives a nominal field of view of 28.8" by 28.8"for the WF Hartmann camera, more than enough for the 20" square FOV ofthe primary detector. This "gross" overfill will shrink appreciably aslonger wavelengths than 500 nm are imaged in the primary beam, i.e. 750nm imaging will be nice and snug with this field of view. Though notexplicitly listed in FIG. 8, the 3.6" spacing plus the plate scale willgive a dimension of 0.45 mm by 0.64 mm on the individual reflectiverectangles on the single checkered mirror 64 and the double checkeredmirror 70, FIG. 6B. Because a one to one finite conjugate field lens 92will be used, the dimensions of these reflective squares will begenerally twice as large on the second stage checkered mirror pieces(these second stage beamsplitters are essentially "splitting" a 4 by 4array of incoming field angles, thus needing a 4 by 4 array rather thanan 8 by 8). These items are not listed on FIG. 8 because they aregenerally derivable entities based on the figures presented in FIG. 8.Turning now to a closer look at the first focal plane 66 and the firststage beamsplitter which will be placed there, a single well correctedfield lens 51 will be placed immediately preceding the first redirectionmirror 62 of FIG. 6A. Its focal length will be chosen to be 500 mm (orwhatever "subtle design variation" conjugates the telescope pupil 34with the first pupil plane 90). Referring now to FIG. 6, the field lenssystem 51 is chosen in the example system to only have a first half 51,and to place an optical flat at the back halves, 53a-d. This is chosenfor many reasons, not least of which to illustrate the "normal"variations which can be made to the "symbolic" design options presentedin the preceding disclosure. Another advantage to this design is thatonly one well corrected lens need be designed and built, not five built.It also should ease optical tolerance. The first field lens (andbeamsplitter, 56) then projects a pupil image onto the first pupil plane90, at which plane a single field lens 92 is placed whose focal lengthin 250 mm (nominal) and whose clear aperture is a minimum of 100 mm (Fnumber<=2.5). Yet another significant variation on the symbolic designof FIG. 7 is chosen, in other words. Rather than having four individuallenses on each of the new beams, normal optical engineering practiceallows the use of one field lens that all four new beams utilize. Theonly real effect of this is to invert the four new beams relative to thesecond focal plane 98 (not a problem). Again, cost and simplicity allowand suggest this choice. The field lens 90 then sets up theaforementioned one to one magnification ratio which can so simplify lensdesigns and mechanical designs. Other choices are quite allowable asusual. The second focal plane 98 therefore has the same plate scale asthe first plane 66, 1"=0.0125 mm. Also, the nominal distance betweenplane 66 and plane 90 is 250 mm, as is the distance from the plane 90 tothe plane 98 (i.e. one to one finite conjugate). At the second focalplane 98 is placed four second stage beamsplitters 100a-d as in FIG. 7.These will also use a single field lens rather than the front and backof FIG. 6. The focal length of this field lens, 51-2nd, is a nominal 64mm placed again directly in front of the first redirecting mirrors ofeach of the beamsplitters 100a-d. This focal length therefore creates asecond pupil plane image a nominal 86 mm beyond the second focal plane98. The lenslet array 52 is therefore placed at this second pupil plane(@52 in FIG. 8). The 86 mm distance between these planes allows thepupil image to grow to a diameter of approximately 8.6 mm at the lensletplate @52. As can be seen, a square lenslet array (both in extent and inthe individual lenses themselves) is chosen with a linear spacingbetween lenses of 240 microns and F numbers of a nominal F/25. Thosepracticed in the art can verify that these figures are consistent withwhat r₀ =7 cm seeing conditions would suggest. The overall size of the"active lenslets" on lenslet array 52 would need to be a nominal 50 mmsquare. A multitude of spot images are then formed onto a 2048 by 2048square CCD with 24 micron square pixel sizes. This CCD was chosenpartially on the basis of choosing a "nominal allocation" by a 5 by 5array of pixels for each and every spot. All in all, this is a bit ofthe chicken and the egg dilemma in choosing these design specs, i.e.what drives what other numbers. In this example, the choice of a 2048 by2048 was rather clear from the original system needs, while the 24micron pixel size (readily available commercially) drove many of theother specifications. The few specifications not mentioned at this pointhave to do with the lateral distances 101 and 103 for both the firststage beamsplitter 56 and the four second stage beamsplitters 100a-d.FIG. 8A has a helpful diagram in this regard. This diagram shows the 16pupil images as formed on the lenslet 52, overlaid on top of the beamdisplacements induced by the beam splitters 56 and 100a-d. As can beseen, an attempt is made to pack these images together in order to use asingle detector 10. This may not turn out to be the end all and be allof packing spots, but it does quite nicely for the example system ofFIG. 8. Referring back to FIG. 8, the lateral displacement of the firstbeam splitter 56 is 14 mm with a nominal 4 mm displacement in the center(slightly larger than the FOV at the focal plane 66), and the lateraldisplacement of the second set of beamsplitters 100a-d is a nominal 7 mmwith the same 4 mm square displacement in the center. It should bestressed that this is a nominal design and engineering considerationsmay dictate small or large changes. In particular, the creation anplacement of the mirror sub-unit(s) 81. FIG. 6 may dictate a largerplate scale and hence a larger center offset. This would ostensibly pushthe pixel area needs beyond a 2048 by 2048 CCD detector. The nominaldesign chosen, however, nicely (very snugly) fills the 2048 by 2048design and shall be used as the explanatory prototype. Hopefully therehave been enough caveats already injected into the discussionsurrounding FIG. 8 to drive home the point that this is not a "strictrecipe" presented in the example system specifications of FIG. 8, but anominal design useful in describing at embodiment of the invention indetails and potentially useful as a springboard into typical engineeringdesign projects.

This concludes the detailed description of the camera system 2, FIG. 1.The disclosure now moves on to a detailed disclosure of thecomputer/data collection system 4, and thereafter into the descriptionof the data processing methodologies.

THE CAMERA CONTROLLING SYSTEM AND DATA COLLECTION SYSTEM, 4, OF FIG. 2

FIG. 9 depicts the camera controlling system and the data collectionsystem. The state of commercial technology is such that this entiredepicted system is readily available by configuring commercialproduct(s). Essentially, it is what is commonly known as a digitalcamera system or a digital imaging system (among many other names aswell). It is well known in the practice of systems engineering (andcomputer system configuration) that there are as many ways to arrangethese commercially available elements as there are systems designers.The figure attempts to do justice to this multitude by leaving therelationships quite general. The only feature which could be consideredsomewhat unique to the needs of the illustrated embodiment of thisinvention are the shutter synchronization sub-system 118 and 40/42. Theoverall system should also be optimized for a high repetition rate ofdigital image acquires, and it would be running two digital cameras, 9and 11; but these are only minor and normal enhancements to a typicalcommercially available digital imaging system. Finally, there is thesoftware 120a-c which oversees and controls the data gathering process.It is customary and typical to utilize commercially available softwarelibrary functions, programs, and compilers, usually sold along with thedigital imaging hardware and host computers themselves, and incorporatethem into the customary system software as auxiliary support programsand user program sub-routines. This disclosure will not get into theparticulars of how these digital imaging systems function, such as howthey direct the processes of reading out the electronic signals from thefocal plane arrays nor how they specifically change analog electronicsignals into digital information and how they specifically store thatinformation, for such details area already known in existing systems.Instead, this disclosure will describe the details which are pertinentto the invention as a whole, and any details which might be somehownovel or specially required by the illustrated embodiment.

For the purposes of specificity, the preferred embodiment of theinvention utilizes CCD detectors as the focal plane arrays of choice forthe visible and near infrared light region. Other cameras operating inthe infrared are certainly possible, as well as non-CCD focal planearrays in the visible/IR. Both CCDs, 8 and 10, are controlled by units108A&B, which in turn are controlled by either a local computer bus in acommercial camera or by a host computer as in a plug-in cameracontroller card (the generic term "Computer bus" 119 will be used andunderstood by those practiced in the art of systems engineering asrepresenting the basic backbone of the data and control pathways). Thereare optional cooling controller units, 110A&B, along with theircorresponding physical cooling units, 112A&B, which are well knowncomponents which can assist in lowering the overall noise levels in adigital imaging system. Here the cooling units can be considered to bewhat are known as common Peltier cooling devices, though other means ofcooling are also utilized. Though the exposure times are typically onthe order of hundredths or tenths of a second, these cooling devices canindeed assist in lowering the effective noise levels of the digitalimaging systems because of the longer readout times commonly associatedwith low noise imaging, and are therefore included in the preferredembodiment. Each CCD, 8 and 10, has its output digitized by analog todigital converters 114A&B and stored in digital frame store memoryplanes 114A&B. The process of digitization typically scans the 2dimensional electronic charge pattern at the faces of the CCDs through ashifting of the electrons which have been generated by the intensitydistribution of the impinging light waves. In other words, an analogsignal comes off of the CCD representing a serial scan, row by row, ofthe entire CCD array. This analog electronic signal is then digitized asper normal electronic apparatus, referred to as an A/D converter, andthe resulting digital values are placed into a 2 dimensional grid ofmemory in an exact sequence at which the electronic signal was read out.For high performance systems, digitizing the electronic signal to atleast a 12 bit precision is recommended [although this is ultimatelydriven by overall system specifications] The frame store memory buffers116A&B are drawn as a stacked array of multiple frames each. This isgenerally a convenience in the assistance of the throughput of digitalimages. where only one frame store buffer per CCD could suffice if needbe. The critical function of the A/D network 114a&b and the framestoring memory planes 116a&b is to generate what is commonly known as adigital image, i.e. an array of numerical values which is a directrepresentation of the intensity distribution of light which existed atthe face of the CCD. The commercial systems which perform this operationare well known to be quite linear devices, that is, the array of numericvalues ultimately residing in the frame store buffers 116A&B do indeedmatch very closely to the relative intensity distribution existing atthe face of the CCD. The last item depicted above the computer bus isthe shutter synchronization and shutter drive circuitry 118, here drawnas a single unit driving both shutters, a driver which opens and closesboth shutters at the same time for the same duration, typically 1/1000thto 1/10th of a second per dual exposure. The overall operation andcontrol of these combined units, 108A&B, 110A&B, 112A&B, 114A&B, 116A&B,and 118, comes from the computer bus (either local or host), andultimately from either software 120a-c or from front panel controls 112,or both. There is also an optional local display device 124, as issometimes common. For the configuration where the camera controller is alocal device hooked to a host. FIG. 9 also depicts a communicationssubsystem 126 and a box representing the host computer 15. It is fair tosay that the combinations of system configurations which could bedepicted below the computer bus line are quite numerous, and it is quitecommon for normal engineering considerations to dictate how to configurethese computer system components. Finally, there is depicted thesignificant option of an image compression and mass storage subsystem,130 and 16 respectively. Depending on many variables, including digitalimage size, number of exposures required per clear image, and others,some systems may decide to immediately compress the recently createddigital images residing in the frame store buffers 116A&B and store themin less expensive slow retrieval mass storage. This is almost anecessity for objects which will require dozens if not hundreds ofexposures without somehow "processing them on the fly" Often, a diskdrive can suffice as the mass storage 16, and as often it could be tapedrives or the up and coming technology of optical storage. Thecompression stage 130 simply allows a greater capacity of images to bestored on a given size storage device without appreciable loss ofinformation. This concludes the detailed description of thecomputer/data collection system 4, FIG. 2. An important end result ofthis description are that this system supports the creation of digitalimages deriving from the detectors 8 and 10, and that there exists acomputing system which can retrieve this digital image information andperform the data processing methodologies of the invention, which is thesubject of the next section.

DATA PROCESSING

To understand the data processing methodology of the preferredembodiment of the invention, it is helpful to first generalize the inputdata to the processing steps. Physically, in the preferred embodiment ofthe invention already described, the input data is what is contained inthe frame store buffers 116A&B and/or the compressed images in the massstorage 16 which once were in these frame store buffers. Additionalinput data is what is commonly referred to as system knowledge: aknowledge of the shape of the primary telescope mirror, a knowledge ofthe spectral characteristics of the light beams 30 and 32, a knowledgeof the pixel sizes on the two CCDs 8 and 10, a knowledge of theeffective sampling window profiles on the individual field angles in thewide field Hartmann camera, and many more items if need be. Referring toFIG. 10, the data input to the processing methodology disclosedhereinafter is:

1) A set 214 of E "distorted" digital images where the distortion isknown to be appreciably variable as a function of position in thedigital image (i.e., field variant).

2) A set 216 of E blocks of data which either directly quantifies(measures) this field variance of the distortion, or else encodes thisinformation in such a manner that subsequent data processing can,through a priori knowledge of the encodation characteristics, derive adirect quantification (estimation) of the field variance of thedistortion.

3) A set 218 of system data, system knowledge, or other environmentalinformation which in any way can clarify ambiguities in the above twoclasses of data or otherwise enhance the overall system performance.Significantly, general statistical knowledge about the distortion wouldbe included here.

The data processing methodology disclosed as the preferred embodiment ofthis invention in effect operates at this most general level ofabstraction.

A processing methodology according to this aspect of the presentinvention creates a single clear image out of a set of E field variantdistorted digital images (or digitized analog images), a set of dataproviding "direct or derivable" information about the field variantdistortion itself, and any other a priori knowledge of the system designor the "statistics" of the distortion itself.

As so defined, this processing methodology is applicable to any linearimaging problem affected by field variant distortions and capable ofsome form of quantifying that distortion and its field variance. Itoperates on a set of distorted images to create a single clear image.This end result is analogous to a machine which would take severalexposed pieces of "blurred negatives" of an object and produce a single"unblurred negative" of that object. An image deblurring system andmethodology, per se. The large dashed boxed in FIG. 10 surrounding theinput 220, the processing engine 222, and the output 224 reinforce thisequivalence to a physical machine operating on film images.

DIGRESSION: GENERAL MATHEMATICAL BACKGROUND

Before describing the details of the flow diagram in FIG. 10, and theensuing further details of the specific preferred embodiment containedin later figures, an extended digression into a useful mathematicallanguage is first in order. This is necessary if for no other reasonthat complicated design considerations can be expressed in concise, butunequivocal terms. A much more thorough discussion of the generalmathematical considerations involved in the basic problem of imagerestoration can be found in chapter 2, "Two-Dimensional Transforms"written by H. C. Andrews, contained in volume 6 of "Topics in AppliedPhysics: Picture Processing and Digital Filtering," ed. T. S. Huang,Springer Verlag 1979 or in chapter 1, "Signal Restoration, FunctionalAnalysis, and Fredholm Integral Equations of the First Kind" written byC. K. Rushforth, contained in "Image recovery: Theory and Application,"ed. Henry Stark, Academic Press, 1987, both of which are available inthe libraries of large technical universities. There are several othersuch sources as well. Rather than continually referencing such sources,I will focus in on those items pertinent to this discovery in thefollowing discussion.

A helpful perspective to view the processing methodology is through whatis commonly known as the linear theory of imaging systems(quasi-monochromatic, unless otherwise indicate). It is summarized inthe boardly accepted mathematical expression:

    g=Hf                                                       (1)

In the general case of this invention where a multiple set of E digitalimages is collected, it is best represented as:

    g.sub.e =H.sub.e f; c=1 to E                               (1A)

Here, f is a vector which represents the unknown object being imaged, gis that set (vector) of digital data collected by a digital imagingsystem, and H is the matrix which describes the linear relationshipbetween each data point in the vector g to each point in the vector f.In actual practice, true linearity is never absolutely achieved, but itis well known that the residual deviation between "linear assumptions"and "non-linear actualities" can easily fall within the scope of a"system error budget."

For so-called "space invariant" imaging systems, i.e. field invariantsystems as this disclosure refers to them, each point in the vector ghas an identical relationship with the set of points in the vector f,except for a shift in indices. Expressed in terms of a digital camera,this implies that one pixel has a certain "view on the object space"while its neighbor has the same view, only positionally shifted [i.e.FIG. 3 where the window function does not change shape as it moves fromone data point to the next.] Many textbooks exist which explore thesebasic points in much greater detail. The point here is that since somany practical imaging systems come reasonably close to this "spaceinvariant" approximation, much of current image processing techniquesrevolve around this assumption since it affords a wide variety ofmathematical and computational simplification. An enormous amount ofsimplification, usually. The deviations of an actual imaging system fromthis approximation are most often small enough so that the errors causedby the approximation are well within system error budgets (the same aswith the "linear approximation" mentioned above). One of the mostbroadly accepted and widely used manners in which this field invariantapproximation is expressed as a mathematical relationship is:

    g=hxf                                                      (2)

or

    G=HF                                                       (3)

These look quite similar to equation 1, but have indeed undergone amajor simplification. Equation 2 is what is known as a convolutionequation (as symbolized by the convolution operator x) and equation 3 isequivalent expression in fourier space, where the capital letterssymbolize the fourier transform of the lower case space functions g,h,and f of equation 2. These are all common expressions in themathematical representation of physical systems and are widelyrecognized and utilized by the engineering community (again, referencethe texts cited above). Even though ultimately ALL imaging systems areFIELD VARIANT, so many systems are close enough to field variance thatthe errors introduced by using equations 2 and 3 in practicalapplications are negligible.

Not so with the problem of imaging objects through a turbulent mediasuch as the atmosphere, however. That is, IF one wishes to image objectswhich are not merely confined to the "isoplanatic region" as previouslydefined. Even when an application calls for the imaging of an objectconfined to the isoplanatic region, it will be found by anyoneexperimenting with the prior art techniques, in relation to thetechniques of this invention, that the "isoplanatic assumption" leads tohigher noise levels (per number of images collected) than thisinvention.

To use the common vernacular of the digital image processing industry,the overall problem of imaging through a time varying field variantdistorting media can be defined as: "Estimate f given a set of E digitalimages, e_(e), each of which has a random (albeit within generally knownstatistical parameters) non-separable space variant point spreadfunction described by the matrix H_(e)."

The very definition of isoplanaticism incorporates the idea that h ofequation 2 does not change other than a shift. As already mentioned, thetransfer function of the atmosphere, i.e. the rows OR columns of H inequation 1, changes quite appreciably across the angles of the sky. Somuch so that even with a shift of angle in the sky of 10 arc seconds ofone degree, in normal "good" seeing conditions commonly referred to asr₀ =10 cm, the corresponding h functions of two angles 10" apart are allbut totally uncorrelated for visible wavelengths. This all buteliminates the massive body of image processing tools which rely on theshift invariant assumption from being employed in comprehensive, globalor wide field solutions to the turbulent imaging problem. This includessimple fourier techniques. In fact, the second and third classes ofprior art systems, adaptive optics and speckle interferometry which aredescribed in the Background section, are both essentially limited bythis field invariant approximation, the first (adaptive optics) due tothe fact that they only measure and correct for paraxial light waves,and the second due to the fact that fourier representations, which ipsofacto impose isoplanaticism, are at the heart of much of the practicedprocessing of their data.

A major foundational novelty of this invention (not here "claimed" perse, but leading to specific methodology inventions) is that the jumpfrom equation 1 to equations 2 & 3 is emphatically not made, and thatother approximations are introduced (which MUST be introduced IF theenormous computational requirements are to be accommodated by less thanall the computational facilities on the planet put together) atappropriate stages, all leading to specific processing methodologieswhich form the heart of the so-called "image deblurring machine" of FIG.10.

Instead, the problem is reformulated into the following:

    g(x)=∫(h.sub.y (x-y)*f(y))dy                          (4)

Equation 4 is nothing more than a restatement of equation 1 (barring thesubtleties of the continuous-discrete transformation process, hereassumed to be "one and the same") with the additionalconsideration--normally quite helpful in actual instrumentation--thatthe coupled variables (x-y) in h are simpler than x and y by themselves.This is to say, the matrix H has a roughly bandlike nature (the word"roughly" is used owing to field variance and geometric distortion) anduse of the coupled variables allows for a simpler conception anddescription of h. The careful observer will note that the "known"statistical properties of a turbulent atmosphere do play a role in thisdecision to utilize the more intuitive coupled variable (x-y), in thatthese statistics indicate a low spatial frequency of the "change inh_(y) " relative to the inherent spatial frequencies in h_(y) ; it is aconceptual "separability" at this stage which will later be transformedinto a quite mathematical one. In our current case of a digital imagingsystem, both x and y can be considered as two dimensional vectorsmetricizing the image and the object spaces respectively. Expanding thisequation out into a rather full representation of a digital imagingcamera, we get:

    g[n]=∫g(x)p.sub.n (x) dx=                             (5)

    ∫(∫h.sub.y (x-y)f(y) dy)p.sub.a (x) dx; n=1 to MN

Equation 5 is the same as that which is depicted in FIG. 3, implicitlyincluding the atmospheric contribution and the telescope optics insidethe field variant h function. Even the pixel sensitivity function, p_(n)(x) is accounted for (again, the monochromatic assumption is in effect).It should also be pointed out that the subscript y on the h functioncould just as well have been an x; the y indicates that use is beingmade of the column vector of the H matrix in equation 1 rather than therow vector. The subscript n runs from 1 to MN, analogous to the discretedigital values that a typical M by N detector would produce. As anaside, equation 5 could be further simplified into:

    g[n]=∫h.sub.n (t[n]-y)f(y) dy; n=1 to MN              (6)

This is nearly a restatement of equation 4 in a discrete format. Thevector t[n] (as a convenience) is a positional vector about which the hfunction can be defined. That is, it preserves the general convolutionalcharacter of the problem. Moreover, it can accommodate lower ordergeometrical and/or "center-of-gravity" "warps" if any advanced designsand methodologies require such, yet retain simpler mathematicalexpressions.

At this point, a kind of axiom is needed which states that much of thesubtlety contained in equation 5, and contained in the use of the vectort[n] of equation 6, is within the proper scope of advanced engineeringand CAN be ignored while describing the illustrated form of theinvention itself. The "errors" which will be introduced into the overallsystem performance will remain just that, errors within an establishedsystem error budget. This is common engineering practice to useapproximations and let the errors fall out where they may. It is thenature of the approximation chosen which will determine whether theseerrors are within acceptable limits. So with this axiom in mind, thefollowing approximation is put forth which can get us back on theexplanatory path of the preferred embodiment of the invention. I returnto a continuous function representation of the problem where the pixelsensitivity function will now be considered to be incorporated insidethe h function:

    g(x)=∫h.sub.7 (x-y)f(y) dy                            (4A)

This is precisely equation 4 only now h tacitly includes a fieldinvariant pixel sensitivity function, p(x). From here, we can follow apath which will wind up with an equation that will inform many of theengineering choices made in the preferred embodiment of the invention.

Focusing our attention first on the continuous function h, one wellknown manner of approximating (arbitrarily closely) h is:

    h.sub.y (x-y)≅Σh.sub.a *Φ.sub.a (y,x-y)(7)

This leads to interesting patterns and insights but not the ones thatthe preferred embodiment is after. Instead, the preferred embodiment canbe better appreciated by employing the following expression which cangenerally be made arbitrarily precise:

    h.sub.y (x-y)=Σz(qs-y)*h.sub.q (x-y)                 (8)

This has been cited in the general literature as a "series expansion ofseparable function," see Andrews cited above. I am currently unaware ofefforts which have gone beyond merely citing the more general form ofequation 8 and noting is lack of usefulness in a global restorationcontext. Here, used with the (qs-y) variable on the z function, carefulexamination will show that this is merely a general interpolationfunction, z, being applied to a "table" of h functions where the tableentry spacing is s. The interpolation function employed is left as anengineering exercise, though the assumption of a bi-linear function is areasonable, practical starting point for this discussion [splinefunctions, standard square box functions, truncated cosine squaredfunctions, etc., are all candidates with various properties] The mainpurpose of equation 8 is the separation of the variables y and (x-y),mixed together in the original h function, into two functions, and in away that retains a simple spatial coordinate system for the fieldvariance (which the grid of table functions with a spacing s provides).The in-operation "typical" frequency spectrums of these two functionswill grossly mimmic the "field variance of h" in z, and the "intrinsiclocal transfer function of h" in h itself. It should be pointed outthat, in the utmost efficiency minded approach to the design of asystem, equation 8 may be somewhat wasteful if the spatial frequenciesof the various basis components of h differ markedly. This is to say,the pole spacing, s, of the interpolants will be dictated by the highestspatial frequency components of h, and may be "overkill" for othercomponents that would have sufficed with a pole spacing longer than s.It is felt that this is of negligible concern at this point and that thematter is entirely left to engineering if the phenomenon truly causesundue inefficiencies. A double summation separable equation 7 could thenbe employed if absolutely necessary. Plugging this new representation ofh into equation 4A, we get:

    g(x)=∫(Σz(qs-y)*h.sub.q (x-y)) f(y)dy=Σ∫h.sub.q (x-y))*(z(qs-y)*f(y))dy                                   (9)

giving a fourier counterpart:

    G=ΣH.sub.q *(Z.sub.q xF)                             (10)

where the shifting, qs, in z is now merely indicated by a subscript q onZ. In general, equations 9 and 10 are no different than equation 1.However, there is now a certain separation of the "localized" (andgenerally lower frequency) components of the field variance from themore familiar transfer characteristics symbolized by the local Hfunction indexed by q. Careful examination of equation 10 will show thatin effect, this is the fourier counterpart to equation 1:

    G=HF                                                       (11)

where the matrix, H, is in italics to differentiate it from in spacedomain counterpart, H. For the sake of computational simplicity andconciseness, it would be helpful to quickly outline the detailedrelationship between the F and G. Expanding F into a series expansionand pulling the summation out front, equation 10 turns into:

    G(k)=ΣF.sub.r *(ΣH.sub.q (k)*(V.sub.q (k)xφ(12)

where it can be seen that everything inside the outer parentheses,including the known and/or fully derivable H_(q) functions, is now asingle known function:

    G(k)=ΣF.sub.r *W.sub.r (k)                           (13)

Once again, this is merely reformulating equations 10 and 11,emphasizing that the rows and columns of H matrix in equation 11 can beeasily calculated on a computer. Detailed examination of equation 10through 13 will find that the H matrix is also roughly banded just asits space domain counterpart H in equation 1. This concludes the general"abstract" mathematical discussion which will be referred to in latersections describing the preferred embodiment in detail.

Turning now to a highly practical mathematical formalism in which tofurther approach the computational necessities of the preferredembodiment, refer back to FIG. 3. The function h in equation 4A isgraphically depicted in FIG. 3, where two particular points in g havetheir own unique h functions projected onto object space. It should benoted that each and every point in g has a unique projection onto(window on) f. FIG. 3 also depicts a variety of steps taken which beginto transform equation 4A into the preferred embodiment of the invention.The main step involves choosing a set of orthonormal basis functionsover which to define f (and subsequently h). It is well known in theengineering community that there are many such sets to choose from, eachhaving differing properties for specific applications, through allserving as a sufficient mathematical description of the functions. Thechoice of basis functions for the preferred embodiment of the inventionis anchored by a uniformly spaced grid whose linear extents along thetwo axes correspond precisely to the dimensions of the pixels on theprimary array 8, with precisely the same ratio of spacing, i.e., thepreferred embodiment uses the most common set of basis functions: thegrid of sinc functions (or delta functions, if we must) which match thedetector. At each "tie point" in the grid--points where the grid linesintersect--there is a two dimensional sinc function (sin(x)/x) centeredover the tie point which is one member of the basis function set. Thecollection of these 2D sinc functions is then the entire set of basisfunctions for f. This is perhaps the most common form of basis functionsused in digital imaging applications (if only implicitly rather thanexplicitly; another common way of describing the very same set ofspatial basis functions is to depict a grid of delta functions in thecontext of a "sampling theorem limited" system, i.e. a Nyquist limitedsystem). Altogether there are MN two dimensional sinc functions,identical to each other except for a positional shift.

Now there are a few practical issues which arise at this point which ameticulous engineer can notice and which, before being swept under therug as prima facie negligible, should at least be highlighted andanalyzed. It is pertinent to the invention, and pertinent to a fullunderstanding of how to implement an embodiment of the invention, tolist these issues and discuss them briefly in the context of building anembodiment according to the present invention. There may arisesituations where this discussion will be highly pertinent to engineeringdetails, especially in systems designed to image through extremelyturbulent conditions.

To begin with, the use of any finite set of basis functions to representsome object 20, FIG. 1, which for all intents and purposes requires aninfinite set of basis functions to be described, brings up a variety ofsubtleties. Physical instruments which attempt to estimate thecoefficients of as many of these basis functions as possible inevitablymust make do with a finite number of data points to work from and hencecannot but form an estimate of the object with a finite set of basisfunctions. Meanwhile, ALL of the object's infinite (orthogonal) basisfunctions in some way affect the data points thus collected, if only ina vanishingly small manner for the vast majority. Most methodologieswhich attempt to form an estimate of the object using a finite set ofbasis functions often experience a phenomena whereby specific basisfunctions "beyond the designed finite number" manifest themselves aspatterns in the data set in exactly the same way as some combination ofbasis functions "within the designed finite number." This is a very wellknown phenomenon in the signal and image processing community, generallyknown as aliasing or as an example of the Nyquist theory (which targetsthe fourier set of sinusoidal basis functions as the working example ofhow aliasing is defined and manifests itself). This phenomenon, and itsgeneral abstract presentation above, is important to this invention notnecessarily at its core, but rather relates to implementation-specificissues. For one thing, this more general presentation can more clearlyshow that errors due to borders and errors due to "frequency aliasing"are two sides of the same coin.

Continuing then, there is the classic and widely utilized theory of thetelescope which culminates in the concept of the diffraction limit. Thistheory, in a rough summary, postulates the existence of a finiteaperture viewing an infinite extent object at an infinite distance insome particular wavelength of the electromagnetic spectrum; Using afourier mathematical representation now, any point P inside an enclosedinstrument will have a window function h (of equation 4) which isabsolutely limited in (angular) frequency to a period of arctan(λ/A),where λ is the wavelength of light and A is the maximum aperturedistance, both measured in the same units. The driving force behind thistheory of the diffraction limit is the idea that, as a luminous pointsource is moved around in the object plane, the complex phase of theincoming planar light waves as projected onto the "pupil" (the open areaof the aperture) can only relatively modulate at this maximum frequencyand no higher. Any point inside the instrument will receive a "nullsignal" from angular frequencies higher than this "diffraction limit"angular frequency with the period as stated above. In broadly acceptedtheory, only the cosine multiplicative component on both the complexamplitude and the complex phase of the incoming light waves createsfrequencies above the diffraction limit, where these two components arevanishingly small at the diffraction limit frequency even for very smallaperture imaging system (they at least, however, provide a solidmathematical negation to the idea of an absolute diffraction limit andto an absolute null space). As a final note on describing conventionaltheory and its telescope design implications, focal plane arrays aregenerally chosen to have a spatial pixel spacing which provides 2 pixelsper cycle of this maximum angular frequency, though some applicationswill specify a greater or lesser spacing for a variety of reasons.Returning now to the relevance of all this to the current invention, itcan be noted that the conventional theory of the diffraction limitposits the existence of planar light waves impinging at the pupil 212.But the entire supposition of this embodiment of the invention, and itsrationale, is that an imaging system will be used in a situation wherethe light waves are not planar but are distorted from a plane wave, andmore importantly, these "distorted plane waves" change as a function offield angle. For field invariant imaging systems, this distortion doesnot really change the concepts embodied in the diffraction limit statedabove, since though the complex phases and amplitudes are distorted,their relative modulation at the pupil plane will be identical asbefore: limited by angular frequency with the period arctan(λ/A). But,once the field variance of the window function h is acknowledged, thesituation summarized in equation 11 once again holds and the plane waveassumption behind the basic idea of the diffraction limit breaks down.Imagine then that incoming light waves from two different directions areseen impinging on an aperture. One just happens to have no distortionwhile the other has distortion. As the field angle is moved along thesky, the phasor diagram of an internal point of the instrument will befound to move at a faster frequency in absolute angular terms than wouldotherwise be the case with the undistorted waves, i.e. at a frequencyabove the "diffraction limit." Furthermore, examination of equation 11indicates that the object spectrum F is being convolved by the(generally symmetric) interpolating function V, thus frequencies "abovethe diffraction limit" contribute to the data vector g.

Some general engineering observations are in order at this point. Theoverall point to the above discussion is that the pixel spacing may needto be somewhat smaller than would otherwise be the case for astraightforward diffraction limited calculation. Furthermore, an objectestimate MAY be sought which exceeds the frequencies of the so-calleddiffraction limit. Having said all of this, it will be found that forlarge telescopes in particular, where the spatial bandwidth of the fieldvariance as reflected in V is much smaller than the spatial bandwidthreflected in H, the effect can be considered negligible. On the otherhand, in smaller aperture telescopes where V begins to appreciablyapproach H's spatial bandwidth, the effect can be significant to thedesign. As a very crude rule of thumb, the pixel spacing should bedecreased from that dictated by normal diffraction considerations by theamount:

    sqrt(p.sup.2 *s.sup.2)/p                                   (14)

where p is the pixel spacing in arcseconds dictated by traditionaldiffraction theory and s is the "table spacing" of the h functions inequation 8. It can be seen that for large aperture telescopes thisfactor is quite negligible in even some of the most severe turbulenceconditions. Be this as it may, the effect is one to take note of,possibly, in non-conventional imaging system design applications.

The above discussion overlaps to no small degree with the issue ofaliasing previously mentioned. The "sale" route is to pick a pixelspacing conservatively above the modified diffraction limit indicated inequation 14. However, this can lead to less photons per pixel andpossibly higher noise levels in certain critical system resultsspecifications. In these cases, often a choice is made to design aspacing which is some significant spacing larger than is dictated byequation 14, improving the S/N ratio at the expense of "tolerable"aliasing. As always, the choice of exact pixel spacing, and thetrade-off between noise and aliasing, is best left to specificengineering refinement.

This concludes the general mathematical description of the basis problemand the discussions on some of the mathematical subtleties that may beencountered.

THE PREFERRED EMBODIMENT OF THE DATA PROCESSING METHODOLOGY

Returning back to FIG. 10, we can now overview this processingmethodology aspect of the preferred form of the invention. We can thenquickly transition to the specifics of the preferred embodiment. Thefollowing discussion is only meant as an explanatory discussion of howthe major pieces fit together, while the recipe for the preferredembodiment will come later.

The data input to the process, 220 as a group, has already beendescribed. The first step in processing the input data is to define aset of "universal" signal to noise ratio segregating basis functions.The engineering choices for exact choice of this set are vast, as anytext on the series representation of functions will attest. This set ofbasis functions is NOT, generally, orthonormal. Indeed, it is generallyan overdetermined set of basis functions, i.e. the total number of basisfunctions is greater than MN (except in the cases of very low turbulenceand in situations where, for whatever reasons, an orthonormal set isrequired; it once again boils down to engineering considerations/errorbudgets). Each member of this set of universal basis functions 226 formsthe basic element where an unambiguous weighted average will take place(I will refer to entities by the number in the box where they aredefined or created). The following step 227 is, in broadest terms, thestep where a user inputs the parameters which will control laterprocessing methods. Here, a choice is explicitly made as to whichreversion methodology will be utilized inside the inner loops 230 and232; the preferred method of the "window function," as will be describedlater, is referenced as the usual choice. The next step 228 in FIG. 10explicitly establishes (allocates memory as it were) both a numeratorand a denominator array for each and every basis function in the set226. These correspond to the usual numerator and denominator of atypical "weighted average" operation depicted generically in FIG. 11.The ultimate purpose of these first two steps is to set up a descriptionof the object estimate, f, which can sort out the signal to noise ratiocomponents of each and every digital image of set 214, and relate thesecomponents between two exposures and between any given exposure and theset of exposures as a whole. Furthermore, as will be explainedhenceforth, use of the weighted average in particular can be used toturn what has commonly been called an "inversion" into what is herecalled a "phase reversion" with post hoc magnitude inversion, or just"reversion." The deliberate phrasing here reflects an appreciation forthe vast work that has been done on "deconvolution" problems, and theimplication of this is that the ubiquitous problem of a possibledivisions by zero or very small values, so common in straightforwarddeconvolution processing methodologies, is all but eliminated by this"reversion" (the remaining subtleties will be discussed). FIG. 11Arelates detailed steps on how a set of E multiple exposures canaccomplish this "reversion" as opposed to the "inversions" which beg forsuch arbitrary solutions. The need for segregating the function f into auniversal set of signal to noise segregating basis functions is alsomade clear by this figure: otherwise there would be an ambiguous andvarying designation of the function which truly is representative of thesignal to noise ratio. Since the reasoning contained in FIG. 11A iscentral to the disclosed method of "undistorting" multiple independentlydistorted images, we turn briefly to describing its steps. Step 1isolates a single member (coefficient) of the universal basis functionset 226 a is the object parameter, b the transform and c the datarealization. Straightforward deconvolution techniques of step 2 lead tothe well known division problem. Step 3 and 4 now include multiplerealizations of the data elements c with the corresponding transforms b,broken out into a phase and magnitude component for reasons which willbe immediately clear. Step 5 points out that s is a perfectly acceptablemetric on the inherent signal level contained in c, assuming the noisestatistics are invariant across the E exposures (if they aren't, justaccount for the additional noise non-uniformity). Step 6 then pulls itall together into a final expression which "puts off" the division untilall exposures have been accounted for. Step 7 illustrates that E, thenumber of exposures, is so chosen as to ensure an average "good" levelof signal by the time this division is made. If a given realizationstill is below a threshold, perhaps it can be dropped entirely. Butbetter yet, a Wiener-type amplification of the magnitudes could be used,as is described in any basic text on signal processing (i.e. thenumerator in 6 could be "amplified" by a signal to noise ratio dependentfactor).

Getting back to FIG. 10, then, the heart of the processing methodologyis the loops indicated by steps 230 and 232 in FIG. 10. The outer loopsteps through all E acquired digital images (and their correspondingdistortion information), while the inner loop steps through each andevery basis function defined by steps 226 and 228. Before the inner loop232 is entered, however, either the input information of data set 216 isdirectly translated into an array named "rev₋₋ array" 235. IF theinformation in set 216 is directly in the form of h equation 4; or anearlier step 234 of transforming "coded" distortion information into theform of equation 4 is performed. The preferred embodiment is of thelatter type, where the input data set is a set of digital images of awhole bunch of spots, hardly in the form of equation 4. Suffice it tosay that steps 234 and 235 are dealt with in great detail in thedescription of the specific preferred embodiment. Also, reference ismade to the "compressed H matrix" in 234 and 235. This language ischosen since it is perhaps the most concise way of describing what isrequired at this stage: an informational efficient encodation of eitherthe rows or the columns of H in equation 1. Because the rows and columnsare variant, BUT relatively low frequency variant, there are hugely moreefficient manners of storing an arbitrarily close estimate of H in thesemore compact manners. Again, in the preferred embodiment of theinvention, step 234 is quite involved, actually, in that it musttransform the spot image data contained in the frame grabber 116B and/ormass storage 16 into the applicable H matrix rows/columns. In this case,the input data 216 is quite "encoded" and step 234 is quite integral.The inner loop 232 is then entered. Step 236 then calculates the arrayof data which will be unique to each basis function member in loop 232;this will be seen to be at the heart of overdetermination, wherein a"smaller" data region has its effects on a much larger region processed.Steps 237 and 238 are then performed which, through choice of theparticular reversion method chosen in step 227, and the consequentassignment of appropriate values to rev₋₋ array 235, calculates thevalue to be added to the accumulating numerator and denominator(elements of 228) corresponding to the given basis function of loop 232(i.e., its corresponding element of the defined set 226). Step 236through 238 are also quite involved within the preferred embodiment. Thevery fact that they are boiled down into a simple accumulatingmultiplication illustrates the close connection between the weightedmean operation and the "reversion" operation. The denominator operationhas the square of the magnitude of rev₋₋ array added into its cumulant,following the principles outlined in FIGS. 11 and 11A; stated in yetanother way, subdata 236 will have an attenuated magnitude approximatelyequal to the magnitude profile of rev₋₋ array, and thus their product inthe numerator necessitates a squared product in the denominator (FIGS.11 and 11A do a better job of describing it, actually). The "next"statements 232B and 230B are then encountered respectively, being theclosing elements of the loops 232 and 230 respectively. Once all digitalimages have been cycled through, anew array called "image₋₋ estimate"240 is created and its initial values set to zero, and then a new loop244 A&B is encountered once again stepping through the entire set ofbasis functions defined in step 226 and established in step 228. Step245 is a generic filtering operation which could be performed at avariety of locations in the processing, even back in step 235; however,this filtering step 245 allows for a chosen optical transfer function towhich the image estimate conforms rather than the "default" equal energytransfer. This is to say, the process without this step would try toboost ALL spatial frequencies to their unattenuated "object space"values, whereas it is generally a better idea to recreate thesefrequencies at the same OTF as, say, the diffraction limitedtelescope--hence the diffraction limited OTF filter would be used atthis step 245. Step 246 then divides the numerator by the denominatorfor each basis function and adds the resultant directly image₋₋estimate. Because the data arrays "subdata" 236 were (or should be, Ishould say) separated in an energy conservative manner, i.e. the sum ofbrightness in all subdata arrays equals the sum of brightness in thewhole image, the additions into "image₋₋ estimate" will ultimatelyremain energy conservative as well. After the loop 244 A&B is finished,the resultant image₋₋ estimate is the "clear image" output 224.

Now that the steps have been generally described, it is easier todescribe one novelty of the processing methodology aspect of the presentinvention. To sum it up in one sentence, the illustrated embodiment ofthe invention simply performs a weighted average inversion--"reversion"as dubbed herein--defined over a generally overdetermined set of signalto noise ratio segregating basis functions. In plain English, theprocess not only "undistorts" the individual images, it does so in amanner that allows for all images to be properly averaged together intoone final crisp image. An individual "undistorted" image IS generallynoisy, and its noise properties non-uniformly distributed. Moreover,this non-uniformity of noise properties changes significantly from oneimage to the next. The present process, through the use of what couldalso be called universal signal to noise ratio preserving basisfunctions, creates a common ground for the image to image noise propertydifferences to be reconciled, while the weighted average "reversion" notonly ensures a maximally efficient signal to noise ratio on the finalimage estimate, it hugely reduces the practical difficulties associatedwith divisions by very small or null entities. Null entity inversion atthe final step 246 would only take place if a given member of the basisfunction set continually finds null values in all its E realizations,and in this case, either the whole member contribution at step 246 canbe thrown out, or a more general Wiener-type filtering can take place atstep 246, where the magnitude "amplification" is related directly to itsinherent signal to noise ratio as reflected in the denominator. [Basisfunction members near the "diffraction limit" generally ensure thatthese latter subtleties are usually present and must be accounted for.]

As a final mathematic note before the detailed recipe for the preferredembodiment is described, we can further analyze why the universal set ofbasis functions defined in step 226 is generally overdetermined. Thisdiscussion leads to further insight on the novelty of this aspect of thepresent invention and its relationship to the prior art. It is wellknown that by employing the standard mathematical technique of thesingular value decomposition (SVD) to equation 1, a solution f can beobtained which is made up of a vector of coefficients multiplying acorresponding series of derived orthonormal basis functions, each ofwhich has a further associated numeric value generally related directlyto its inherent signal to noise ratio. Rather than get into the detailsof this process, we only need note here the "non-uniform noiseproperties" alluded to above can find precise mathematical expression inthe SVD formalism without any need for "approximation." Every student ofthe SVD processing methodology knows that the matrix representing theorthonormal basis functions and the vector representing the inherentsignal to noise ratio both change according to the H matrix alone, notits data realizations. It should not be surprising, then, to discoverthat for a "typical" random media changing the H matrix as it does, itis quite difficult or downright impossible to find a single"independent, orthonormal" set of basis functions which can serve as theperfect segregator of signal to noise properties across multipleinstances of the H matrix. The transformation from any given instance ofthe SVD solution onto this hypothetical "universal" set of basisfunctions will force a scrambling of the otherwise "clean" signal tonoise vector (sigma, in the SVD vernacular). It is important to choose aset of basis functions, in step 226 FIG. 10, which can appreciate theimplicit or explicit "approximation" which is taking place between the"SVD ideal" and the chosen universal basis function set--that is, thisshould be done such that a given threshold of acceptable error ismaintained. Computational efficiency and resource usage efficiencyalmost demands that approximations be made, while the overall systemerror budget provides the limit to these approximations. As alluded toearlier, for very low turbulence situations and for system designs withliberal error budges, a good candidate for an independent orthonormalset of basis functions (i.e. normal and maximally efficient) is aspatial weighting grid of lower spectral resolution sinusoids, e.g.:

    Φ.sub.a,b (x)=(sin(cx)/cx)*sin-&-cos (dx)              (15)

where c^(z) c^(z) d^(z) d=MN (this is cryptic, admittedly, but the wordsshould speak for themselves and it is quite immaterial to the disclosedinvention). It should also be noted that in the case where the overallimaging system can be well approximated by the isoplanatic assumption,i.e. equations 2 & 3 can be profitably used, then it can easily be seenthat the sinusoidal basis functions of the standard fourier componentscan fully serve as the "universal" set of signal to noise segregatingbasis functions. Indeed, this very concept is superfluous in that case;field variance is what gives it real meaning. It should also be notedthat the prior art, most notably in the paper "Deconvolution fromwave-front sensing: a new technique for compensating turbulence-degradedimages," J. Primot et. al., J. Opt. Soc. Am. A Vol 7 No 9, p 1598, notonly implicitly makes use of this clean segregation of the signal tonoise components in the isoplanatic case, it use of the time average ofthe correlation by the time average of the auto-correlation, equation 3,is a degenerate case of the general weighted average reversion listedabove, i.e. it happens to degenerate into global orthonormal basisfunctions which are derivable through less general operations such ascorrelations. It is clear, however, that the correlation andautocorrelation operations lose their efficacies when applied to fieldvariant systems and thus do not extrapolate beyond the "isoplanatic"assumption. In any event, this functional congruence of this equation 3to the more general principles outlined in FIGS. 11 and 11A dictate thatthis invention disclosure will not presume to claim the use of thisspecific processing technique within "fully isoplanatic" methodology. Onthe other hand, "fully isoplanatic" applications which utilize "fieldvariant technique enhancements" certainly fall under the purview of thisinvention and will be claimed accordingly.

THE DETAILED PREFERRED EMBODIMENT

The description of the preferred embodiment picks up where it left offwith the detailed descriptions of the camera system 2 and the computercontroller 4. This preferred embodiment also makes use of the quitespecific formats and variables of FIGS. 1 through 9. However, it shouldbe clear to those practiced in the art that based on the above fewsections, a much wider variety of input/output conditions could beaccommodated by the processing methodology of the invention.

FIG. 12 contains a schematic of the information which is best describedas belonging to the input system data 218 of FIG. 10. Exact numbers areused in order to illustrate their relationships and in order to specifya typical system which can be built using current technologies. FIG. 12depicts a simplified schematic of the pixels 300, represented as points,of the digital images derived from the digital camera 9, with p" 302being the equivalent field angle spacing between pixels (0.02" as perFIG. 8, the example design). A square image format is chosen with atotal number of pixels N 304 on each of the axes, where FIG. 8 has N at1024. The xs 306, implicitly overlaid directly on top of certain pixels(or fractions of pixels as the case may be), represent the nominalcenter points at which the pupil function has been measured by the widefield wavefront camera 11 (i.e. the pupil function for a given fieldangle). These nominal center points 306 have a nominal field spacing oft" leading to a pixel spacing T 308, where in the example design of FIG.8, 3.6" leads to T=180. The light reaching the primary camera focalplane array 8 and detected by it has the spectral characteristics (asshaped by the overall system, the CCD 8, and the dichroics 28 and 58)contained in the inset FIG. 12A, where the central wavelength is 500nanometers and the spectral curve has a fully width half maximum of 25nanometers (in this example; remember that this FWHM figure isnegotiable relative to overall system error budgets). FIG. 12A also hasan example of what the spectral sensitivity curves on the wavefront beammight look like both for a broad band system, the full curve, and for abeacon optimized system (see later sections) as the dashed curve.Harkening back to the discussion on the extension f the spatialfrequency coverage beyond the diffraction limit, as symbolized inequation 14, we discover that the ratio of extension is trivially smallin this design and so the nominal pixel spacing of p"=0.02" will indeedbe chosen. It should be noticed again that the example system of FIG. 8will have its pupil function sampling extending beyond the extent of thedata pixels; this is a convenience as already mentioned and need notnecessarily be the case. Instead, this field angle sampling could befully inside the extent of the data pixels and use would be made ofsimple extrapolation methodologies for the pixels at the outer edge ofthe detector. The choice of a field angle spacing at 3.6" should bequite sufficient to measure turbulence conditions at and somewhat"worse" than the well known seeing parameter r₀ =10 cm, probably down tor₀ =7 cm or even lower with appropriately more error incurred. Aspreviously discussed, there is no absolute wall at which point finalresults fall off a cliff in an effort or signal to noise sense. A fieldspacing of 3.6" would even give reasonable results down to a seeinglevel of r₀ =5 cm, though various error and noise mechanisms wouldcertainly begin to show up in the final results by this point. Any moredetailed discussion on this point belongs in the theory sections and innormal engineering refinement.

Referring back to FIG. 10, the first bona fide processing step is the"non-computational" definition of the universal S/N segregating basisfunction set. FIG. 13 schematically describes the set chosen for thespecific example described in FIG. 8. First, the digital image is splitup--using square spatial weight functions--into a spatial grid of "datasub-zones" 320. Here, the weight function used is a 32 by 32 pixelregion of unity value with zeros everywhere else, stepped across theentire digital image space. At the expense of computation time andresources, this could as well have been a bi-linear weighting functionor a smoother interpolating function instead of a 32 by 32 squareregion, but the example of the square region suffices for many overallsystem error budgets. It can be noticed that now the digital image datais split up into 32 by 32 regions, each with 32 by 32 pixels. This issomewhat a coincidence: i.e. the 32 by 32 within a 32 by 32; since a 2Kby 2K detector would give a 64 by 64 region of 32 by 32 pixels. It alsoreflects a preference for powers of 2 which will become more apparentshortly. At this point, each new data region, of which one in particular322 is highlighted in FIG. 13, defines its own unique projection ontothe solution space of f, here defined as a 256 by 256 point grid with acorresponding 256 by 256 array of fourier components (its own UV plane).The position of this 256 by 256 spatial plane relative to the solutionspace of f is not shown, but it can be seen as shifted inside thatsolution space so that its center coincides with the center of the datasubzone from which it derives. It will be the fourier basis functionswith the coefficients defined in the UV plane 328 which will be utilizedas the components of the universal basis functions defined in step 226,FIG. 10. Since each of the data regions 320 defines its own unique basisfunctions, it can be seen that there will be a total of 32*32*256*256,or 67108864 individual basis functions (64 Meg), or a factor of 64larger than would otherwise be needed to fully specify the solution f.This is to say, this set of basis functions is overdetermined by afactor of 64. Several engineering points must necessarily be insertedhere in order to explain this.

First of all, there are clearly choices which could be made which willgreatly reduce the overdetermination factor. The example presented inFIG. 13 chooses a conservative approach more appropriate to empiricallyderived error budgets as opposed to well-defined theoretical ones. Thethinking behind the choice of a 32 by 32 region is that it is "wellwithin" the "isoplanatic" conditions such that the pixels within theregion can be treated in a normal isoplanatic processing methodology.Likewise, it is not too small as to unduly increase the computationaltime and resources in the hunt for diminishing returns (negligibledecreases in error and noise). It is clear that the ideal universal setof basis functions would be to have each pixel define its own 1K by 1KUV plane, giving 10¹² individual basis functions and anoverdetermination factor of 10⁶. The choice made in FIG. 13 is thereforea compromise which, significantly, includes the desire to keep theprocedure within a simpler framework not only to explain it more easily,but to simplify the real world software programming requirements whichso often includes non-mathematical experts, or at least persons who arenot experts in all technologies involved in the illustratedimplementation. It has already been pointed out that weight functionsother than a square 32 by 32 region could be used to split up the dataregion, also. It is quite conceivable, but to my present knowledgebeyond practicality, that a basis function extension of the individualdata regions could be found which is maximally efficient ANDanalytically smooth in describing the "reversion" process, but: will itbe computationally efficient as well?. These are questions for future"efficiency refinement" R&D. On a rather abstract plane, the path ofengineering refinement can be summarized as "minimize the computationaltime/computational resources for a given level of error budget";mathematically, find: MIN (DIM(#BASIS @ TURBULENCE), EXECUTION TIME) @ERROR₋₋ LEVEL; The DIM alludes to a minimum dimensional entity,correlated directly to memory needs (# of basis functions for a giventurbulence) and partially correlated () to execution time, which canfully represent the reversion procedure for a given error level; but itmust be chosen so as not to overly burden the execution time above thehighly efficient FFT procedures of the preferred embodiment, (which willbe described). This sort of exercise is a classic engineering refinementtask.

Returning now to FIG. 13, both the choice of a 32 by 32 sub-zone, andthe choice of a 256 by 256 "halo" 324, were largely determined by theknowledge of the (very gross) statistical properties of the turbulencein the atmosphere; this statistical "knowledge" being generallycontained in the input information 218 of FIG. 10. The choice of the 32by 32 pixel data sub-zone 322 was already shown to approximate a regionwhere the model of isoplanatism could hold, i.e. where equation 2 and 3could profitably be employed within the overall system error budget (IFeq 2 & 3 are needed; for more refined cases they are not needed, e.g.,in reversion methodologies which work within local low order varianceconditions). Furthermore, this can lead to the use of efficientprocessing methodologies such as the two dimensional fast fouriertransform, to name the hands down favorite. The choice of the 256 by 256"halo" 324 could just as well have been the entire 1024 by 1024 solutionregion in extent; the truncation to a smaller region (based on a powerof 2) merely limits it to a region of significance, that is, thetruncation of the area largely reduces the computational load whileintroducing a "negligible" error. The definition of "negligible" fitswithin basic engineering considerations, and this example decided thatgoing with a 128 by 128 "halo" would not have been negligible (forturbulence levels approaching r₀ =7 cm and worse). A general rule ofthumb that can be used in choosing the size of the "halo" 324 is basedon the well known "long exposure" profile of a point source. This ideais quite well known in the industry and has long been used to define theinherent resolving power of unassisted telescope viewing through theatmosphere. It is a common rule of thumb that in "seeing conditions ofr₀ =10 cm, using the 500 nanometer wavelength, a full width half maximum(FWHM) of 1" is generally found for the "long exposure" profile. Thelong term exposure profile represents the statistical weighting functionthat any given data point has on object space. A choice must be made to"ensquare" some percentage of energy of that profile. If we look at thefield angle sampling spacing of 3.6" (308, FIG. 12) and the 500 nmcentral wavelength, it is clear that the overall system should workreasonably well down to seeing conditions of r₀ =7 cm and perhaps a bitlower (higher turbulence levels) even. In this case, we can expect anFWHM of 1.5" to 2.0" at worst, or 75 to 100 pixels in width as thecurrent numbers have it. The 256 by 256 square halo, on the other hand,translates to a conservative 5.12" by 5.12" square region on the objectspace. Subtracting the 32 by 32 inner breadth of the data region, thisgives a 112 minimum pixel extent from an outer pixel of the 32 by 32region to the edge of the halo, depicted as item 326 in FIG. 13. Theidea here, verified by mathematical/empirical models and common sense,is that any given pixel "can be affected by" the object space containedin this halo, and therefore any "solution projection" from this smalldata region would affect the object estimate to at least this extent.The areas of the object space beyond this halo (the x% of light energybeyond this) will then show up as error signal, or "aliased" signal,within the overall system error budget. The choice of 256 by 256 in thisexample ensures that this "aliased signal" is in the few percent regionand lower, down to the nominal design point of r₀ =7 cm. The reason forthe earlier "abstract" discussion on the concept of aliasing was toaccommodate the idea that "frequency aliasing" and "border aliasing" aretwo intuitive sides of the same coin.

To finish the description of step 226, FIG. 10, which is furtherdetailed by all of FIG. 13, we can not that the 32 by 32 data sub-zone322 is centered inside the 256 by 256 halo 324. The corresponding UVplane 328 is there drawn as a 256 by 256 complex array (if onlyimplicitly). It is more useful at this juncture to be more explicitabout how to regard these entities as computational elements. Inparticular, FIG. 13A contains the basic in situ encodation of theseuniversal basis functions as they CAN be used on a computer (as ispreferred by myself and as used in the preferred embodiment). Here,instead of envisioning the 32 by 32 data sub-zone 322 as centered in thehalo, it is instead seen as 4 (16 by 16) data zones at the four cornersof the 256 by 256 spatial halo. This is certainly not necessary but doeshelp to simplify the UV plane representation by placing "the center ofthe data sub-zone" over the origin of the implicit fourier transformspatial domain. There are alternative ways to do this and many referencematerials to choose from in helping to decide which one to use. Alsodepicted in FIG. 13A is the explicit recognition that the spatial domainis a real, as opposed to a complex, function and therefore only requireshalf of the UV plane to be fully represented, the other half merelybeing the mirrored complex conjugate (as is well known in the industry)of the displayed half. Finally, this process step keeps track of whichdata sub-zone is which (330), therefore being an indication of therelative shift that needs to be applied to the individual basisfunctions in comparison to the global basis functions of the solutionspace.

Turning now to FIG. 14, we find the full and specific process steps--orflowchart--for the software which carries out the standard operationsdepicted in FIG. 10. As with the case of the specific embodiment of thewide field Hartmann design in FIG. 8 in relation to the general designin FIG. 7, the flowchart in FIG. 14 will make specific choices ofoptions outlined generally in FIG. 10. Prototype software and generalimplementation strategies have made use of the common software languageknown as "C," and therefore much of the presentation of the flowchartsconform to the basic concepts employed within that programming language.In particular, the use of callable "sub-routines" are used which thencall forth other figures and descriptions, as opposed to describingthese sub-routines in absolute linear execution format. The use ofsub-routines also tremendously eases the discussions on the variants ofsub-routines, since the discussion on variants can be confined to thoseappropriate sections and not overly burden the description of the mainprogram itself.

The preferred embodiment comprises a computer program written, debuggedand compiled which shall be called "main" 350, as is quite archetypicalto C language computer programs. It is well known within the professionof software engineering that such a program as this could easily be asub-routine inside a larger program, such as a graphical to text-baseduser interface program or a larger project organization program. Apartfrom whatever normal software engineering considerations which enterinto the creation of computer programs, the first programming step 352explicitly listed here is the generic step of declaring variables,opening input and output files, and otherwise allowing for the inputfrom the user of variables which will modify the ensuing behavior of theprogram. The next step 354 is to allocate and clear (apply zeroes asinitial values) the three arrays which correspond to the complexnumerator and the real denominator of the universal basis functions 328FIG. 13A. The preferred embodiment specifies these as four dimensionalarrays of 2 byte signed integers, though the DENOM array could just aswell be 2 byte unsigned integers since it will only contain positivevalues. These arrays could also be 4 byte floating point numbers or even8 byte double precision values, where here they will be used as "encodedintegers" meaning that implicit or explicit scaling values will bemaintained relating the stored digital values to initial grey valuescales in the input image data [This is a common software techniquewhich eliminates the need to store values as floats merely forconvenience; this is to say, there will be little to no information lossin storing the values in as compact a form as possible. In the casehere, it can be seen that the array sizes can get quite large (i.e. 192Megabytes), and a storage savings of 2X can be helpful].

These three arrays can be interpreted as a 32 by 32 array of 128 by 256UV plane elements: that is, the real numerator 356, the "imaginary"numerator 358, and the magnitude denominator 360. Implicit in thesearrays is the fact that the Nyquist frequency has been thrown out, theopinion being that it is not worth the special programming attention italways seems to demand (including many "x+1's" on array indices). Themeanings and formats of these arrays have already been described. Thenumbers indicated show that these are three 64 Megabyte arrays giving atotal of 192 Megabytes. Relatively speaking, this is a large amount ofmemory needed and it is quite common to set up these arrays, in whole orin part, on "longer access" memory than the usual "fast main memory" ofa computer. The choice has little relevancy, other than normalengineering and cost considerations, to the rest of the program. Thenext step is the beginning of a classic for-next loop 362. This loopcycles through the "E" digital sets of "speckle exposures" provided asinput 214 FIG. 10. The implied character of such a loop is that when thelast element of the loop is reached at the "next" statement 362B, theprogram resumes directly after that next statement [as is well known inthe industry] Bearing in mind the nature of how a loop performs,therefore, the first programming step inside the loop 362 is a globalwavefront data reduction sub-routine 364 which will be describedfollowing the description of the main program. This global wavefrontsub-routine 364 will be seen to have a choice of several methodologieswhich in general trade off accuracy with computation time (and otherthings as well, certainly). It is a normal function of engineering tooptimize this trade-off. In general terms, the function of thissub-routine is threefold. First, it sorts out and processes the "digitalimages of the spots" 216, FIG. 10, into an array of complex pupilfunction estimates corresponding to each and every nominal field anglesample 306 FIG. 12. Second, it computes the "point spread function"corresponding to the field angle center points of each of the datasub-zones, using both pupil function "interpolation," the spectralcharacteristics of the lightwave 30 FIG. 4 as additional information(again, this will be described in detail later). Finally, it furtherprocesses this array of pupil function estimates into forms which thenext loop 366 can with utmost efficiency utilize. The details of thislast function of step 362 will be discussed in the description of thewavefront sub-routines. The output of step 364 is therefore either a setof digital values directly relating to an array of point spreadfunctions indexed by position as in FIG. 12 (equivalently, its UV planerepresentation), or a set of digital values which further processesthese numbers in manners to be described later.

The next step is what is commonly known as a nested loop or an innerloop 366. Here, all 32 by 32 sub-zones (320, FIG. 13) are loopedthrough. It is a matter of software engineering preference whether thisis viewed as, AND programmed as, a 32 by 32 array or as a 1024 linearvector. Depicting this loop with only one statement 366 implies thattreating it as a one dimensional vector can simplify FIG. 14 itself. Thenext step 368 is to clear a 256 by 256 array of, either, 2 byte integersor a 4 byte floating point values, the choice being dictated by storageconsiderations, execution time (for conversions mainly), andavailability of appropriate "in-house" of "commercial" software for theensuring 2D FFT routines. This 256 by 256 array in step 368 will bereferred to as "subzone₋₋ array" 370. The next step 370 places thedigital values contained in the current (as defined by loop 366) dataregion 320 FIG. 13 into either the center of subzone₋₋ array 370 as inFIG. 13, or into the four corners of subzone array as depicted in FIG.13A. The choice is purely a matter of personal preference with someamount of housekeeping necessary in either case. FIG. 15A schematicallydepicts the operation of step 370, with the choice of loading the datainto the center of subzone₋₋ array chosen for graphic simplicity. A twoby two array of numbers rather than the full 32 by 32 array is also usedfor graphic simplicity. The next step is a conventional 2D FFT performedon subzone₋₋ array, where a very wide array of commercially availablesoftware can be utilized (all with varying software formats and dataformats, but all within the normal scope of software engineering). Theoutput (denoted by the ≧ symbol) of this operation will either directlybe a 128 by 256 array of complex elements 274 (two digital values) orcan be transformed using normal mathematical processing into such. Here,this new array is "conceptualized" as one array called "uv₋₋ array" 374,but as a programming element, it is stored as two arrays: uv₋₋ array₋₋real[128][256] and uv₋₋ array₋₋ imag[128][256] At worst, a commercialpackage will require a second 256 by 256 array of zero valued elementsto be the "imaginary components" of the input image. An FFT output isthen given as two 256 by 256 arrays, the real and the imaginary uv_(')array 374 can then be constructed by considering only half of each ofthe two arrays, i.e. a 128 by 256 half of both the real and theimaginary array. Again, the "Nyquist" component is thrown out; or, itcan be stored in the imaginary component of the DC element. All of thisodd subtlety is quite well known in the art of image processing. Theonly practical point is to keep track of what one's own preferences areand what the 2D FFT routines require.

Next, step 376 calls the local wavefront sub-routine 376, whose outputis a 128 by 256 array of complex values here called "rev₋₋ array" 378,also breaking out into rev₋₋ array₋₋ real[128][256] and rev₋₋ array₋₋imag[128][256] rev₋₋ array has precisely the same format as UV₋₋ array,only it will have different individual digital values determined in thelocal sub-routine 376. As noted in the discussion on the globalwavefront sub-routine 364, the local sub-routine 376 also has a varietyof methodologies to choose from which are (generally) intimately matchedto the choice made for the global sub-routine 364. As such, both ofthese sub-routines will be discussed immediately after the discussion onthe main program in FIG. 14. For our current purposes, it is sufficientto understand that the output array, "rev₋₋ array" 378, is a UV-planerepresentation of the complex phase "reversion" that will be applied tothe local data contained in sub-zone array 370. VERY GENERALLY SPEAKING(to be clarified in the discussion on the wavefront sub-routines, andwhich was already clarified in the earlier sections on the mathematicalbackground), rev₋₋ array contains the information necessary to reversethe phase distortion in the data, while also carrying information aboutthe magnitude distortion which will be used in the ensuing weightedaverage [essentially the information on each UV plane element'sintrinsic signal to noise level properties AND its magnitude"attenuation"] For the purposes of the main program, however, we needonly realize that rev₋₋ array and uv₋₋ array are identical in format andframe of reference.

A third "nested loop" 380 steps through each and every (two-valuedcomplex) element of the matched arrays UV₋₋ array and rev₋₋ array. Inother words, it loops a total of 128×256 times, or 32768 times. Thefirst step 382 inside this loop 380 is to perform a complexmultiplication of the current complex value in uv₋₋ array 370 by thecorresponding complex element in rev₋₋ array 378, where current isdefined by the index on loop 380 and the complex multiplication isobviously performed on the real and imaginary components of thosearrays. The resulting complex value can be designated as the matchedpair (RE, IM) 384, being floats, doubles, or integers as the case maybe. Step 386 assigns the variable "mag" 388 to be the sum of the squaresof the two values RE and IM, i.e. mag=RE² +IM². It should be emphasizedthat this is the square of the complex magnitude of the complex variable(RE, IM) 384. Step 390 converts whatever working data format (ints,floats, etc.) was used during the preceding calculations into valuescompatible with the three arrays NREAL, NIMAG, and DENOM; 356, 358, and360 respectively. Step then adds RE to the existing value of NREAL[ . .. ], IM to NIMAG[ . . . ], and MAG to DENOM[ . . . ] It is understoodthat the index values on the three arrays correspond to the currentvalues defined by the loops 366 and 380. The exact operation in step 392corresponds to the "C" language "+=" operation, i.e. NREAL[][][][]+=RE;etc.

After step 392, we encounter a cadence of "next" steps for the loops380B, 366B, then 362B. After the last iteration of loop 362 iscompleted, i.e. after all digital exposures have now been cycledthrough, the pertinent result is that the three arrays NREAL 356, NIMAG358, and DENOM 360, are now completely filled with the numericinformation corresponding to the complex numerator and real denominatorof the weighted average on each and every universal basis function 328FIG. 13. The next step 394 is to allocate and clear a new array to benamed image₋₋ estimate[1024][1024], here explicitly being a twodimensional image estimate (spatial domain) with array sizes matchingthe primary detector array size. Again, clearing means placing zeros asinitial values for all elements of the array. This clearing is importantsince there will be an ensuing accumulation process expecting nullvalues to begin with.

The next step is to enter a loop 396 which once again steps through the32 by 32 array of data sub-zones. Indeed, loop 366 could have been usedhere with an "if last exposure, then" statement in front of the "next"statement 366B, but this current arrangement is more clear. The nextstep 398 is a new nested loop, stepping through each and every UV planeelement of the current data subzone. It can be seen that loop 396 and398 merely step through all indices on the NREAL, NIMAG, and DENOMmatrices. Step 400, inside the nested loop 398, once again makes use ofthe (temporary) 128 by 256 complex array "uv₋₋ array" 374, assigning thereal element, uv₋₋ array₋₋ real, with the division of NREAL[current] byDENOM[current] and assigning the imaginary element, uv₋₋ array₋₋ imag,with the division of NIMAG[current] by DENOM[current] There may or maynot need to be data type casting involved if the division is betterperformed by floating point arithmetic and NREAL, NIMAG, and DENOM arenot already in that format. Also of note in step 400 is the addition ofan element called "filter[][]" which is essentially a scalar quantityreferred to in box 245 FIG. 10. In cases where a Weiner-like filteringof the final image is preferred, the division by DENOM becomes purely afunction of DENOM as is the common Wiener method in the image processingcommunity. The closing "next" statement 398B of loop 398 is thenencountered. The result of completing loop 398 is that uv₋₋ array 374now contains the UV plane representation of the weighted average of thecurrent data subzone's basis functions. Effectively, this is the current(indexed by loop 396) data subzone's contribution to the final imageestimate 224, FIG. 10. [As an aside, the image estimate f could just aseasily have been a 512 by 1024 complex array, i.e. its UV planerepresentation. In that case, step 404 about to be described would havebeen a complex multiply rather than a spatial shift. As usual, thespatial domain approach was chosen not only for simplicity, but also tobetter understand the effects of data borders.] There is one practicalproblem with uv₋₋ array 374 at this point: detailed analysis will showthat it represents a "full" deconvolution, meaning that it attempts tobring ALL frequency components up to Nyquist (or the diffraction limit,if you prefer) back to their original magnitude. In other words, it istrying to estimate the object up to the absolute informational limitsset down by the solution grid spacing and image extent. What is muchmore appropriate at this stage, however, is to merely attempt toestimate the object to the resolution level that the underlyingdiffraction limited telescope would permit as if there were no obscuringatmosphere. This means that uv₋₋ array 374 should be "filtered" using a(digital) Airy disk (or the annular equivalent; or whatever).Ultimately, the choice of filtering is wide open. Filtering uv₋₋ array374 with the natural optical transfer function of the underlyingtelescope is merely the clearest nominal choice. It could as well be afully symmetric gaussian distribution with a spatial FWHM in the sameneighborhood as the FWHM of the telescope's inherent Airy disk.Therefore, step 401 performs this filtering via a multiplication in theUV plane, as is quite common. The multiplication is generally with areal symmetric function such as a gaussian of the OTF of the underlyingtelescope. Importantly, the ultimate purpose of filtering at this pointrather than fully at step 410, is that otherwise there will be a muchhigher signal content "wrapping" around the 256 by 256 boundaries, i.e.,the signal level does not subside enough by the time it reaches theboundaries in the case of the unfiltered uv₋₋ array 374. This can createadditional error in the final image estimate f 394. Step 402 thenperforms an inverse FFT on uv₋₋ array 374, turning it into a 256 by 256spatial domain representation "subzone₋₋ array," 370, where we make useof that earlier (temporary) array name. The inverse 2D FFT f routineusually would be a part of the same off the shelf type library routinesused in step 372. Loop 404 then adds subzone₋₋ array 370 into image₋₋estimate[][], taking care to add it while referenced to the center ofthe data subzone, throwing out values which exceed the image estimateborders, AND taking care to unshift it out of the form depicted in FIG.13A, IF that was the form of storage chosen. [The format of FIG. 13A mayseem somewhat awkward at this point, but discussion on advancedmethodologies in creating "rev₁₃ array" 378 sheds some light on how thisshifted form can be useful.] FIG. 15B graphically depicts loop 404. Inthis figure, the 256 by 256 array "subzone₋₋ array" 370 can be seencentered over the original data region from which the subzone itselfderives, hence the offset 408. This all is within the context of theentire 1024 by 1024 image₋₋ estimate. An additional considerationderiving from the presence of a data boundary is also alluded to in thepicture and deserves a brief tangent.

It can be seen in FIG. 15B, and noted by step 406 inside the loop 404,that the portions of the shifted subzone₋₋ array which fall outside ofthe 1024 by 1024 boundary of the image estimate are merely ignored. Thisdisclosure refers to this manner of dealing with data boundaries as the"fundamental spatial domain boundary method." It is truly a point ofdeparture for a wide range of philosophies and methods which can bebrought to bear on this well known problem. As a point of comparison,had I chosen the fourier representation of f and the complex multiply inthe UV plane as may "shifting method," I would have wound up with whatthis disclosure refers to as the "fundamental frequency domain boundarymethod," wherein the parts of subzone₋₋ array which are "thrown away" inthe spatial domain method would be "wrapped into" image₋₋ estimate in amirrorlike fashion on opposite sides of the 1024 by 1024 boundary ofimage₋₋ estimate. In a crude sense of conservation of energy, thefundamental fourier method is preferable, but the spatial domainapproach was chosen here once again for its ability to explain the basicboundary problem, its implications, and to provide an intuitive base toderive preferable (though generally more complicated) optimizedsolutions. To state the basic problem to intuitive English, objectswhich are inside the estimate boundary of f will have some fraction oftheir light energy fall outside of the physical detector 8, FIG. 1, witha generally larger fraction as an object is closer to the projectedgeometric boundary of that detector; the information carried in thislight energy is lost. Likewise, objects outside the projected geometricboundary will have some fraction of their light fall within the physicaldetector 8 and will, in general, generate eventual "aliased" signalcomponents (using quasi-informational terms). The implications of thisproblem can be conveniently split into two categories closely aligned tothe underlying applications of the imagery itself. We can refer to theseas 1) quantitative and 2) qualitative. In the approach toward creating"better" solutions to the boundary problem than the "fundamental" one,there are some common denominator elements and some unique elementsrelative to these two basic categories. We can scratch the surface ofthese for the sake of completeness of this disclosure, but the bulk ofrefinement once again falls under normal engineering (and scientific, tobe sure) refinement. The simplest and most intuitively applicabledescription on the basic implications of the boundary problem isgraphically contained in FIG. 16. This figure is deliberately"non-quantified" other than the pixel extent numbers. It merely has aplot of "general error" as a function of linear distance of a solutionpixel from the solution border. A note is made on how corners will havecompound effects. At a distance of 128 pixels, half of the extent of ouruniversal basis functions, the error reaches a uniform level solelyrelated to the fixed finite size of the universal basis function.Crudely speaking, the curve is related to the integral of the longexposure point spread function. At least from the qualitativestandpoint, a rule of thumb "implication" falls out from this wholediscussion on the boundary problem: Consider the solution region withina distance of the long exposure FWHM (full-width half maximum) from theborder to be corrupt and throw it out. Thus, in the example of thespecific case preferred embodiment, the FWHM is 1 arc-second in r₀ =10cm seeing, or 50 pixels in extent on our 2.5 meter telescope system at500 nm; so we would consider the "qualitatively acceptable" solution tobe 924 by 924 pixels in extent (50 pixels subtracted on all sides,corners subtleties ignored). The second drawing in FIG. 16 depicts thisgraphically. Individuals experienced in the imaging art will see thisfor the "arbitrary" definition that it is, but at the same time, mayaccept the value of having simple rules of thumb such as these in orderto simplify the bigger picture. One final note is that a practicalproblem emerges of the difficulty that now exists in maintaining a powerof two format for the data. Subsequent data analysis will now have 1024by 1024 of "valid" data to work with. It is quite a tempting option,then, to have a primary detector which is slightly larger than 1024 by1024, say 1280 by 1280, and then always shoot for a "valid" data regionof 1024 by 1024 (for instance). This then shifts the power of two"problem" to earlier stages of processing in FIG. 14, AND brings up thedifficulty of finding commercially available detectors, frame grabbingRAM, etc. [The "power of two" desire is driven both by software andhardware considerations.] Fortunately, this is all well within the realmof normal (but tedious) engineering and science. The disclosure of theinvention will leave the topic of the boundary problem at this muchsaid, the essentials having been presented.

Finishing up with the main progress and FIG. 14, the loop 404 and itsmain processing component 406 adds the values of subzone₋₋ array intoimage₋₋ estimate, with the "offset" 408 determined by the current indexon loop 396. Step 404B closes the loop, thereupon encountering the nextstatement 396B for the loop 396 which is cycling through the datasub-zones. The entirety of loop 396 can be seen as stepping through eachdata sub-zone, finding its weighted average contribution to the finalestimate, then adding that contribution to an accumulating image₋₋estimate 394. After the loop 396 is completed, there is an optional(additional to 401) filtering step. Using this "post-filter" is quitecommon, and in the context of the preferred embodiment of thisinvention, can be helpful in attenuating anomalous signal caused by theartificial truncation of the data field into 32 by 32 pixel regions,likewise with the use of a 256 by 256 basis function extent, and evenpossibly residual effects of the "optical" sampling artifacts caused bythe fixed wavefront sampling by the Hartmann camera. These "errorsources" once again should be placed into the context of an overallsystem error budget, thereby driving decision IF and HOW to use thisoptional post-filtering step 410. These error sources can be"attenuated" by other means obviously, such as through refined choicesof universal basis functions, refined methods of computing rev₋₋ array378, and so on. We finally reach the output, then, a "diffractionlimited" image if we must use the vernacular. Again, practitioners ofthe art will realize that the use of the term "diffraction limitedoutput" is arbitrary, and other more technical ways of describing theoutput would be: "highly improved strehl ration image" or "greatlyenhanced OTF imagery" or "quassi-Wiener filtered image" (if S/N ratioshelp determine the ultimate transfer profile), and so on. The last itemmentioned here, the "Wiener" aesthetic, is the one that I generally useat least at the conceptual foundations, it is the one that relates thecore of the results. Nevertheless, the broader astronomical community ismore comfortable with the "diffraction limited" terminology or someequivalent.

THE WAVEFRONT ESTIMATION SUB-ROUTINES

There are two reasons that these sub-routines are split out into theirown section. First, they are generally quite specific to a Hartmann typewavefront camera in particular and to automatically distorted imagery ingeneral (as opposed to the most basic category of "distorted digitalimagery"). Second, they are relatively involved and contain their ownalternative methodologies (generally based on the typical "accuracy vs.resources" engineering dynamic) and therefore they would have dilutedthe discussion on the main program to the point of unintelligibility.Now that the rationale for this has been stated, we can turn once againto the specific preferred embodiment which is based on Hartmann "spot"data being input data 216 in FIG. 10. Recall, also, that our specificdesign of FIG. 8 posits on 8 by 8 array of sampled field angles. Onefinal note is that those experienced in the art of "traditional"paraxial type Hartmann wavefront sensors are aware of a wide range ofprocessing methodologies and other practical techniques that aregenerally available in the public literature. There is no attempt madein this disclosure to go into detailed discussion on these methods.Instead, reference will be made to items which belong to this wealth ofpre-existing art, and only those items which are unique to this form ofthe invention, or those items which are part of the prior art but whichhave an intimate relationship to those items which are unique to thisform of the invention, will be discussed at length and in detail.Overall reference is here made to Southwell, "Wavefront estimation fromwave-front slope measurements," J. Opt. Soc. Val. 70 No. 8, p998-1006and to Herrmann, "Least₋₋ squares wave front errors of minimum norm," J.Opt. Soc. Am. Vol 70 No. 1, P28-35 as excellent treatments of thebasics, where each of these articles in turn have several otherpertinent references. Concerning more recent "refinements" to the basictechniques involved, reference is made to Cho & Peterson, "Optimalobservation apertures for a turbulence-distorted wave front," J. Opt.Soc. Am. A. vol 6, no 11, p1767-1775.

FIG. 17 depicts a simple overview of the wavefront distortionmethodology. It is specific to the use of the Hartmann camera asmentioned. It is important to point out that the end product of thisentire wavefront section and the described sub-routines is to create thecomplex array "rev₋₋ array" 378, FIG. 14. The wavefront methodology canbe split into three procedures as depicted in FIG. 17, the first ofwhich being further split into MacroCalibration and MicroCalibration asrubrics. The calibration operations are essentially independent from themain program outlined in FIG. 14. It is a procedure which (ideally) willbe performed only as often as system performance requires: perhaps oncea day, or perhaps once a month even for MacroCalibration, but certainlyafter particular physical changes have taken place with the overallcamera system 2, FIG. 4. The calibration procedure as a whole outputs afile of calibration information which will be sued by the two remainingprocedures, both of which are sub-routines in the main program of FIG.14. The overall purpose of the calibration procedure is to determine thenominal (Macro) and estimated instantaneous (Micro) centers (nullpositions) of all the Hartmann spots, find out which field angle in thearray of field angles that each spot belongs to, match up each spot toan nominal location in the pupil plane, and to determine at what pointis a given spot "on or over the practical edge" of the pupil plane andnot getting enough signal level relative to the other spots. This lasttask is referred to as finding the "good spots." The second procedure ofthe three in FIG. 17 is call the global wavefront sub-routine and hasalready been discussed in general within its role as step 364 of themain program 350, FIG. 14. FIG. 17 restates in overall purpose and theensuing discussion will fill in the details. Likewise, the thirdprocedure listed in FIG. 17 corresponds to step 376 in FIG. 14, andlikewise, it will be discussed in further detail. As can be seen in FIG.17, the three procedures comprising the overall wavefront distortionmethodologies are arranged in increasing frequency of usage in thenormal operation of the illustrated embodiment. The distinction madebetween these three already crosses the line into engineering efficiencydetail, that is, it splits out the overall wavefront sensingrequirements into initial stages of computation time management. Weshall now discuss each procedure in turn.

Of all the three procedures listed in FIG. 17, the calibration procedureis most affected by which type of Hartmann camera system is beingemployed. The splitting of the calibration into two delegates the bulkof the work to the Macro routine, while the Micro routine "fine tunes"the estimates of the instantaneous spot center references. Threeparticular choices of Hartmann camera system design were detailedearlier in FIGS. 4, 6, and 7, being the simple configuration wide fieldmodel, the 2X model, and the "full" model respectively. Conceptuallyspeaking, the task given to the Macro Calibration procedure in FIG. 17is simple. It can be stated as: For each and every spot contained in an"typical" digital image from the Hartmann Camera system (input 216, FIG.10):

1) find which angle it belongs to, and to which nominal area in thepupil plane does it correspond.

2) find it nominal "averaged" center location (or zero-slope reference)as a function of the pixel grid on CCD 10, where the center is found tofractions of a pixel spacing--as will be discussed.

3) determine is its nominal center point on the pupil plane is tooclose, or over, the edge of the pupil support, thereby precluding a spotfrom being considered a "good spot"; that is to say, its inherentsensitivity will fall below some pre-determined threshold.

Additionally, each field angle in the array of sampled field anglesneeds to be matched to particular location on the primary CCD 8, FIG. 4.We can refer to this broadly as the positional correlation of the twoinput data sets 214 and 216 of FIG. 10.

Using these three tasks as guidelines, and the need for positionalcorrelation as the parallel companion, much of the macro calibrationbelongs to the realm of common sense and basic engineering. There isobviously a great deal of regularity in the arrangement of the spots.Contrarily, there may also be subtle deviations from this regularitywhich, especially in the case of the full Hartmann camera system of FIG.7, would require a rather involved software scheme for "keeping track ofthe spots." It should be pointed out that prior art (paraxial) Hartmannsensors already must deal with, and have already come up with finesolutions for, a somewhat simpler form of macro calibration. It is"somewhat simpler" due to the fact that there is only one field angleinvolved, where in the wide field case there are interlaced fieldangles--multiply interlaced in the 2X and full designs. Referring backto FIG. 4B, it contains the simplified version of the interlaced fieldangles for a 4 by 4 field angle sampling Hartmann camera system (simpleconfiguration). The idea of spots "fading out" as opposed to abruptlystopping is roughly drawn in as well, suggesting the need to thresholdthe "good spots" from the bad ones. [Again, paraxial prior art alreadyhas developed solutions for this type of problem. Also, yet anotherreminder should be made that the clear circular fade out in FIG. 4Bwould actually be more of a fade out of the squares as a unit, not theindividual spots within the 4 by 4 sub-groups.]

The preceding discussion and FIG. 4B are quite sufficient to map outspecific software embodiments. Be this as it may, a description offurther considerations is in order, and detailed description of an"inefficient but thorough" output file will be presented in order tomesh sufficiently with the discussion on the micro calibration and thetwo remaining procedures depicted in FIG. 17.

It might be possible to design the Hartmann camera system with suchtight opto-mechanical tolerances that much of the macro calibrationcould be performed "by design." This is to say, the physical tolerancewould unambiguously assign spots to field angles and pupil positions,and field angles to primary CCD (8) pixel locations. This wouldgenerally involve micron level tolerances for some components, such asthe CCD relative to the lenslet array, and slightly eased tolerance onother components. It would be yet more difficult to accomplish this forthe 2X and the full models. It is equally probable that this kind ofultra-precision--for this particular "need" at least--is not necessary.What is more needed is in-operation rigidity rather than absoluteaccuracy. Even this need for "rigidity" is contingent on whether or not"engineering methods" such as the micro calibration are available andeconomical. In any event, one extremely simple method for assigningspots to field angles (as a first task) is to "scan" a bright star inand around the paraxial field angle of the telescope. Those experiencedin the art of telescope operation will understand that, assuming wegather nominal digital exposures, the bright star will variously "lightup" the interlaced arrays of spots belonging to the individual fieldangle in which the star currently is positioned. Provided that thebright star is positioned such that it is sufficiently centered in agiven field angle, those spots belonging to that field angle readilymanifest themselves. If it is not centered, either no spots light up, ortwo "sets" of spots will light up such as in the full model of the widefield Hartmann design, when the star falls at the border of two (ourfour!) regions. Once a bright star unambiguously lights up a singlearray of interlaced spots, standard prior art methodologies can treatthose spots as if they were from a normal paraxial Hartmann sensor, findtheir nominal centers (over an ensemble of exposures), find their"fade-out" characteristics, and assign an offset to the spots relativeto the known pupil plane support. It can be seen that this lastprocedure accomplishes all three tasks listed above. Again, once thespots are identified which belong to the same field angle, prior artHartmann techniques can be used almost verbatim. It remains to establishpositional correlation between the field angles and the primary CCD.

During the above "scanning of a bright star," the primary CCD 8 shouldbe recording a "centroid record" of the position of the bright starduring each exposure in the series of exposures taken during the"scanning." Realizing, then, that what we are now looking for is thenominal center point of a given field angle as projected onto theprimary CCD, we realize that there are a variety of ways we canestablish this estimate. Perhaps the simplest (conceptually, at least)method would be to isolate a given field angle and to make very smallsteps to position of this bright star within this field angle. Two datapoints are collected during each exposure: the position of the centroidof the bright star on the primary CCD, and the total brightness of allof the spots belonging to that field angle. [The centroid is part of thebread and butter of the prior art in paraxial Hartmann sensors. Briefly,it establishes the center of gravity of the brightness distribution of agiven array of digital values, where that array of digital values isusually peaked and quasi-symmetric. It typically finds a resultant"center" to fractions of a pixel, i.e. one tenth of a pixel or, as acurrent practical limit, one hundredth of a pixel.] After randomlymoving the star in order to create a sufficient ensemble of data pairsin this fashion, a new quasi-symmetric profile manifests, and any numberof methods can be used to determine the center of this distributionincluding another centroid-like operation or a "peak-finding" operation.The end result is that each field angle [n,m] (n=1 to 8, m=1 to 8 in ourspecific 8 by 8 example) is assigned a sub-pixel resolution position onthe primary CCD. It should be clear to those practiced in the art that awide range of engineering refinements can be applied to the aboveprocedures, most definitely including manners of software efficiency andmanners of deducing overall patterns rather than doggedly measuring eachof the 8 by 8 sampled field angles.

Let us now clearly define the required output of the macro calibrationroutine, realizing that the above discussion is sufficient to outlinehow this output is generated. The us of the "C" struct variable type isenormously helpful in this context. The output can be seen in its mostthorough (but preferably inefficient) form as:

OUTPUT CALIBRATION FILE 420

struct field₋₋ angle [8][8] {

float x₋₋ CCD₋₋ position, y₋₋ CCD₋₋ position;

int number₋₋ spots;

float x₋₋ Hartmann[number₋₋ spots],y₋₋ Hartmann[number₋₋ spots];

float x₋₋ Pupil[number₋₋ spots],y₋₋ Pupil[number₋₋ spots];

}

The reason that this is "inefficient" is that no allowance is made forthe regularities involved, regularities which could be exploited forreduction in storage and computation time. Nevertheless, this disclosureshall use this as the foundational approach. The output is arranged asan 8 by 8 array of structures, each with the associated elements insidethe brackets. Each field angle has a nominal center point on the primaryCCD stored in the first floating point values listed. The second linecontains one integer value which gives how many "good spots" were foundfor this particular field angle. Note that this allows for each fieldangle to have a slightly different number of good spots than theensemble average; i.e., field angle X may have 120 good spots, whilefield angle Y has 122 good spots. The third line contains the x and ycoordinates of the nominal "ensemble average" center location of all thegood spots for this field angle. It has "num₋₋ spots" number of floatingpoint values referenced to integer pixel locations on the CCD 10 FIG. 4.The issue of finding the best "center" of these spots is clearly amatter for normal engineering; what we can say about it here is that,ultimately, these center points will represent the "zero wavefrontslope" condition which is well known in the prior art of Hartmannsensors. As with these prior art sensors, there will always be need fora reference from which to determine the instantaneous wavefront slope.The use of ensembles of spot image realization across many field anglesincreases the statistical resolution of the ensembling technique, inrelation to system which use reference beacons to help track theinstantaneous null positions of the spots. The last line in the outputstructure is the corresponding nominal center of each spot in theconjugated pupil plane. This too is an x/y floating point array ofelements, "num₋₋ spots" in size. The indices on these arrays and onthose of the preceding line correspond to the same spots. The referencescale of the pupil plane is rather arbitrary, where one can easily usestandard metric units and assign the x/y axes in any orientation thatone pleases. The only practical matter involved with scaling is to keeptrack of the relationship between the pupil reference scale and theprimary pixel reference scale, relative to the primary beam centerwavelength; but this matters little at this point.

FIG. 18 takes a closer look at the situation obtaining for the specificpreferred embodiment. It can be noticed that the top drawing in FIG. 18corresponds to the appearance of a digital image deriving from the fullor "4X" beam splitting Hartmann camera of FIG. 7. In particular, itdepicts the (nominally optimized) configuration as would appear for the2.5 meter telescope, utilizing a 2048 by 2048 CCD array as the specificfocal plane array 10, FIG. 4 (i.e. the example system, see FIG. 8 oncemore). Calculations will show that there will be a little over 1000"good spots" per field angle, giving a total of over 4000 floating pointvalues per field angle in the above calibration structure, with thestructure itself being an 8 by 8 array. This winds up being about 1Megabyte of calibration information if our floating point entities are 4bytes in size. [A reasonable size; besides, remember that goodengineering practices can hugely reduce this by exploiting theregularities.] The middle drawing 430 blows up a small section from oneof the pupil-images-turned-to-spots from one of the 16 beams. The bottomdrawing 432 then further blows up this middle drawing, focusing in onmean iso-intensity lines for two neighboring spots overlaid on a pixelgrid pattern (using the nominal 5 by 5 pixels per spot allocated in thenominal design). Note that these intensity distributions arepolychromatic according to the spectral properties of the Hartmann lightbeam and the CCD sensitivity, FIG. 12A. They are also intended torepresent non-turbulence high brightness conditions, whereas FIG. 4Balready showed the in-operation movement of the centers of gravity ofthe spot images. Their "FWHM" isolines (50% @436) are seen to be a bitwider than the pixel sizes and their 2%-of-peak isolines appreciablyoverlap. This is drawn arbitrarily here please remember, where eachspecific design may find different isobrightness results; i.e., theseare as usual drawn for expository purpose. The purpose of highlightingthe overlap of the brightness distributions is to illustrate that "spotcrosstalk" will in general be a more central issue in the wide fieldHartmann cameras than in the design of prior art paraxial Hartmannwavefront sensors. In fact, it is the spot crosstalk issue which drivesmany of the practical limitations of the wide field Hartmann camera ingeneral. In particular, it limits the fineness with which the fieldangles can be sampled. The fact that the FWHM is comparable or largerthan the pixel size is quite open to engineering refinement as well.FIG. 18A, very much like FIG. 4B, then depicts a hypothetical inoperation distributions of the spots in a high brightness situation,except this figure draws the iso-brightness lines 436 [as opposed toFIG. 4B and the centers of gravity. These lines are referenced to thepeak brightness at the center of the spot]. The vectors drawn are fromthe quiescent nominal centers (the +'s 434) to the instantaneous centerof the spot. This vector has been generally accepted as being linearlyrelated to the wavefront slope existing at the lenslet subaperture (seeSouthwell previously cited), and hence becomes the core informationcarrier of the phase-difference measurements of the complex pupilfunctions (exactly as with the prior art paraxial systems). Thisin-operation situation, and the issues of spot cross-talk and low photoncounts, will be discussed in the two remaining sub-routine discussions.It is the purpose and methodology of the micro calibration which can nowbe discussed. Very simply, the micro calibration "fine tunes" theestimation of the quiescent center points of the spots relative to themore coarse estimates created during the macro-calibration. The purposeof this micro calibration and its methodology are essentially equivalentto prior art paraxial Hartmann sensors. As such, only an overview willbe stated here. After a telescope system has been configured accordingto FIG. 2, and after the macro-calibration has been performed, theinstrument will of course be utilized on real objects. The movement ofthe telescope from one object to another, the "small" relative movementsof the components of FIG. 2 as an object is "tracked," and the change ofinterchangeable components giving rise to component movement, allcontribute to small changes in the position of the nominal non-turbulent(zero wavefront slope) spot centers from those found duringmacro-calibration. It is thoroughly impractical to re-performmacro-calibration for these small changes, thus, some method of"macro-calibration" may or may not need to be performed, all dependingon how these movements play out in terms of the system error budget. Therange of solutions that can be brought to bear on this "need" formicro-calibration is again wide, with much to be outright borrowed fromparaxial Hartmann sensors systems. This disclosure will present yetanother simple minded methodology, leaving more detailed solutions forengineering studies. The idea behind this micro-calibration is that thefull wide field Hartmann design of FIG. 7 can be considered sufficientlyrigid for most applications, and thus the only movement of nominal spotcenter locations will be a global "x-y shift," or offset, which will beapplicable to ALL spots. Thus, those practiced in the art will recognizethat mis-calibration will only affect the global tip-tilt estimate ofthe pupil function (the second and third zernike coefficient). If agiven Hartmann system is found to have more than this simple form of"error" between the nominal spot centers and the true instantaneous(quiescent) centers, then these can be modelled as zernike coefficientand treated the same as the tip-tilt is about to be treated. Assuming,then, that it is only the global tip-tilt correction which is sought inthe micro-calibration, it would be possible to keep a global average ofthe spot location across a large ensemble of short exposures, say 100exposures of more. The macro calibration would have given its own globalensemble average. As the system operates, this "floating average" isalways compared with the macro-calibration average, and the differenceis always accounted for as the "micro-calibration" to the nominal spotcenter locations as given by the macro-calibration stage. The ensembleshould be large enough so as to not add its own tip-tilt noise to theongoing system performance, but small enough to meaningfully measurewhatever slow drift there is in the tip-tilt parameter. Another moresophisticated method of micro-calibration would use a bright isolatedstar somewhere in the field of view as a kind of anchor point in alarger feedback loop in the main program in FIG. 14. However, thesekinds of solutions go beyond the scope of this disclosure and arelargely connected with optimizing overall system performance, ratherthan enabling operational performance. It can be seen that both of theseabove ideas would run in a parallel fashion to the main program in FIG.14, feeding updated "offset information" to all of the spot centervariables contained in the (third line) of the above calibration datastructures. This disclosure will leave the micro-calibration issue atthis much said, suggesting that the literature on paraxial Hartmannsystems be consulted for further ideas.

We can now turn to discussing the second procedure, "step 2," in FIG.17, the "global wavefront sub-routine" 364 of FIG. 14. During thediscussion on the main program 350, FIG. 14, it was stated that thisglobal wavefront sub-routine has three general tasks to perform.Reiterating, they are: 1) sorting out and processing the "digital imagesof the spots" 216, FIG. 10, into an array of complex pupil functionestimates corresponding to each and every nominal field angle sample 306FIG. 12; 2) computing the "point spread function" (and/or its UV planeconjugate) corresponding to each of the center points of the datasub-zones 320, FIG. 13, by "interpolating" the pupil functions onto thefield angle corresponding to the center of each data sub-zone, alsousing the spectral characteristics of the lightwave 30 FIG. 4, depictedin FIG. 12A, as additional information; and 3) further processing thisarray of point spread function estimates into forms which the loop 366,FIG. 14 can with utmost efficiency utilize.

FIG. 19 can help better explain how these general requirements areimplemented. The first step 500 of FIG. 19 is a non-processing stepwhich clearly describes the input data to the overall global wavefrontsub-routine. Recall that this sub-routine is called for each and everyexposure in the main program, FIG. 14. Thus the input data can be seenas 1) a single "digital image of spots" such as is depicted in the topdrawing of FIG. 18, and 2) the calibration data structure 420 asdescribed earlier. [From here on, it will be assumed that if amicro-calibration procedure is utilized, it will function in parallel tothe main program and will effectively update the calibration datastructure as required, and the global wavefront sub-routine need onlyreference the stored center references and be assured they are thelatest estimates available.] The first processing step in FIG. 19 is aloop 502 through the 8 by 8 array of sampled field angles. The ultimatefunction of this loop is to determine the complex pupil functionestimate at the nominal center points of each and every field angle inthe 8 by 8 array. [Advanced approaches to the problem would not pose thesituation quite in these terms; instead, they would pose the task as"deriving a new set of spatially filtered but discrete" complex pupilfunction estimates. The distinction derives from the tension betweenease of explanation and setting up system optimization. This disclosurewill continue to err toward the former (ease of explanation) since thereis no loss of "enablement" by doing so, and the advanced optimizationscan merely be seen as post hoc refinements, described after the basicenablement is finished.] Step 504 does pay some lip service to thepotential parallel micro-calibration, where whatever updates arerequired to the zero wavefront slope spot center references, the nullpositions, (third line is calibration data structure 420) can be appliedat this point. Step 506 begins a nested loop through the number ofspots, "num₋₋ spots," for the current field angle as given in thecalibration structure. The first step 508 inside the nested loop 506 iseffectively a small sub-routine, but we can discuss it as one operationhere since it is derived from common sense. Essentially, it finds theintegral pixel value of the closest "peak" to the nominal center pointof the current good spot (remembering that the nominal center is definedin the calibration data file structure). The "peak" can be defined as "abrightness value over a certain threshold with its neighboring pixelbrightness values all lower than itself" (as I said, common sense;moreover, prior art Hartmann sensors routinely utilize such processes asthese). Those practiced in the art of image processing will realize thatthere are certainly subtleties involved with defining this peak andconditions which would confound certain algorithms which search forlocal peaks. Complete spot overlap of neighboring spots, identicalbrightness values next to each other, AND local maximum caused by lowphoton count noisy data are common conditions which come to mind.Nevertheless, the preferred embodiment of the invention merely requiresan algorithm that will "find the nearest peak" whether it is the "rightone" or the "wrong one." The engineering idea is simply "design the peakfinding algorithm, and more importantly, the nominal spot separationdimensions (of the physical Hartmann design), such that X % of the timethe peak finding algorithm finds the correct peak for a given level ofturbulence." X % should typically be in the 99% region or above. For ourpurposes at this stage in the disclosure, step 508 merely finds theclosest peak and calls it good--where all the overall system performancerealizes that it will be wrong a small percentage of the time: (100-X)%. Step 510 then performs a "local centroid" about the peak pixel thusfound in step 508. A local centroid processing steps are then describedin FIG. 19A, where a nine pixel neighboring is utilized in thecalculation (the peak pixel and its 8 adjacent neighbor pixels). Thevariable "low" 474 refers to the lowest brightness value among the ninepixels. The variables px[] and py[] in step 482 can take on the values-1, 0, or 1, according to their x and y positions about the centralpixel. The resultants, xcenter and ycenter, will be found withinfractional pixel position precision, and are returned in reference tothe offset of the central pixel of the input 3 by 3 array. Thosepracticed in the art will once again notice that spot cross-talk cansignificantly affect the results obtained by this operation for thatpercentage of spots whose "instantaneous centers" come within a twopixels of each other. Indeed, the cross-talk effects are everpresent.This disclosure will have an extended exploration of this topicimmediately following the current section on the core wavefrontsub-routines, FIG. 17. For now, we will accept whatever "error" isinduced by step 510, FIG. 19, and the 9-point local centroid as part ofthe system error budget. Step 512 then loads the two (temporary array)variables grad₋₋ x[] and grad₋₋ y[]with the respective differences ofthe measured centroid and the nominal "zero slope" center. The temporaryarrays grad₋₋ x[] and grad₋₋ y[] both have "num₋₋ spots" number ofelements corresponding to all the good spots of the current field angle,and the current refers to the loop 506. Step 514 is an optional systemdiagnostic step, where a warning flag "warning1" is incremented (itshould be cleared before loop 502 is entered) is the root of the addedsquares of grad₋₋ x and grad₋₋ y is larger than, say 1.5 pixels in the 5by 5 pixel/spot nominal design case. Later, the value of this warningflag is an rough indication of how turbulent the atmosphere is and howmuch spot cross-talk error might be occurring. An alternative is merelyto record the mean square deviation for all spots; whereas the techniquestated here is more indicative of explicit spot mis-identification. Step516 is optional (though recommended) and is common in prior art paraxialHartmann sensor systems. The basic thinking behind this correction isthat the object giving rise to the spot data may itself have abrightness distribution whose "center of gravity" does not coincide withthe exact center of the nominal field angle being sampled. If we imaginea bright star 1" away from the nominal field angle being sampled, itwill nevertheless be contained in the square steradian field of view ofthe field angle sampling window (of the 4X design) and will skew thecentroid accordingly. For the specific preferred embodiment presented inFIG. 8, it can be seen that this "bright star" can bias the centroid by1/3rd of a pixel, a very appreciable amount when we consider thatcentroid precisions of 1/20th of a pixel are routinely achievable. Thus,step 516 presents an optional method of estimating the "C.O.G." (centerof gravity, to use the common descriptor) of the object within thecurrent steradian field of view of the nominal current field angle, andusing this information to "correct" the grad₋₋ x[] and grad₋₋ y[] valuesfound in step 512. Those practiced in the art will see this as another"global tip-tilt" correction to the given complex pupil function of thecurrent field angle; and indeed, those practiced in the art probablyhave been using some form of this correction on prior art Hartmannsystems. The added importance of this correction to the wide fieldHartmann design device from the fact that the tip-tilt error becomes averitable field angle dependent "wobble" in later system performance,rather than a somewhat arbitrary global shift as it is in the paraxialHartmann designs. Though the advanced methods of wavefront sub-routinespresented in FIGS. 26A&B, 27 and 28 will quite explicity deal with this"problem," we can sketch in a simple solution even at this juncture.Specifically, the ensemble of speckle digital images in the input dataset 214, FIG. 10, can be added together and the resulting image filteredby a suitable estimate of the long term seeing point spread function(e.g. a 1" FWHM quasi-gaussian). A given field angle will be multipliedthis image by its own average steradian sensitivity profile (e.g. asmoothly rounded 3.6" by 3.6" square region in the specific preferredembodiment). The centroid is found for this newly created data set, andthe difference between it and the nominal center of the field angle isaccounted for via step 516, FIG. 19. There is no point in being moredetailed about this particular procedure since the advanced methodslater discussed have more elegant ways of dealing with this potentialerror source; again, it is also common in prior art Hartmann systems.The "next" statement for loop 506 is then encountered in FIG. 19.Putting aside the question of the "error sources" contained in thearrays grad₋₋ x[] and grad₋₋ y[] which will be explored later on, theresultant of this nested loop 506 is identical to what would obtain froma prior art Hartmann sensor system. It should be no surprise then thatthe next step encountered, 518, simply labelled as a prior art"wavefront estimation sub-routine," adding the "least squares" as aclearly identifiable trait. That is, step 518 turns the phase differencemeasurements (combined with knowledge of the positions of those phasedifference sample points relative to the known aperture supportfunction) into an absolute phase estimate with an arbitrary constantphase term. Reference is directed in particular to Southwell's 1980paper for the mechanics and considerations of this sub-routine. Thosepracticed in the art can borrow (probably verbatim) the current routinesthey are already using for "phase difference to absolute phase"conversion programs. As usual, the constant phase term is neglected,i.e. the integration constant is immaterial. It is fair to say that thebasic principles of this phase estimation process are well known at thispoint and that current work involves evolutionary refinement (see paperby Cho and Peterson, for example). One particular point that does havepertinence at this point is that much of the earlier literature (but byno means all) posits the need to perform this process 518 "very fast" inorder to keep up with the physical mirror actuators as in prior artadaptive optics systems. It is important to realize that this need for"very fast algorithms" is emphatically not a requirement for the basicinvention disclosed herein--it will be a consideration, however, for"realtime" implementations of this invention--and this fact will beexploited in the advanced methods in the section on engineeringoptimization. Step 519 is an optional wavefront estimate refinement stepwhich will be extensively described in the sections surrounding FIG. 27.It is written as optional since indeed it is, but the preferredembodiment of the invention will typically perform this function. Forthe purposes of discussing FIG. 19, this step can be summarized asattempting to mitigate spot cross-talk effects as well as attempts toutilize higher order zernike coefficient estimation on the spots if suchinformation would be useful in better estimating the phase[][] Step 520is yet another optional refinement step, a step which generally has beendetermined to be "negligible" in most formulations of the basic problemof turbulence correcting imaging system. This step is an optionalsub-routine which performs a least squares fit to what has been called,among many other things, the scintillation of the pupil function; whichmerely means that the whole complex phase function (magnitude and phasetogether) are estimated rather than merely the phase alone as in step518. Though I don't have explicit studies to cite, part of thebackground folklore of adaptive optics systems in particular is thataccounting for the scintillation function removes a degree of errorwhich is dwarfed by practical error sources deriving from physicalinstrument noise sources and by errors resulting from residual "phase"after step 518 is completed. Be this as it may, this optional step 520is presented here for several simple reasons. The main reason is that noone can say if and when other error sources may "practically" subside tothe point where the scintillation function error would become anappreciable performance enhancement. Another simple reason is that itnever hurts to empirically experiment even when solid theoreticalreasoning suggests you are wasting your time; you never know what otherinteresting artifacts may pop out, such as cross-correlation statisticsbetween phase and amplitude information to the point that thescintillation information can help refine the phase information (longshot? I won't presume to guess). In any event, it is conceivable that asub-routine place at step 520 could "integrate" the photon count at eachspot, and a (generally low(er) order model) scintillation function couldbe least squares fit to the resulting grid of "intensity" values. Spotcross-talk and spot mis-identification are clear difficulties to dealwith, as well as polychromatic effects. This disclosure will leave it atthis much said and leave whatever "discoveries" there are to be made forfuture R&D AND on-site empiricism; this disclosure will also touch uponthe scintillation issue in the engineering refinement section.

Whether step 520 is performed as a bona fide scintillation estimator, orwhether it merely changes the incoming single element phase array into atwo element complex phase array with constant magnitude 1.0, theresultant output after step 520 can be seen as a (variable offset, onefield angle to the next) grid of complex values representing the complexpupil function at the current field angle defined by loop 502: herenamed "complex₋₋ phase[][]," implicity having a real element and acomplex element, generally floating point, and defined on the pupilsupport function such as annular support of the specific design examplein FIG. 2. There are several things we can now state about this array.These are: 1) It is a flawed estimate of the "true" complex phasefunction existing at the precise center of the nominal field angle,owing to a) it is spatially filtered (as a function of field angle), andb) it has "noise"; and 2) It exists on a "grid," or more generally on aset of basis functions, which are unique to this particular field angledue to its unique relationship with the lenslet aperture; that is, ithas a variable offset because of the 16 beams which give rise to the 16circles of spot images depicted in the top drawing of FIG. 18, eacharray of spots has a quasi-random relationship with the lenslet array(i.e. an offset from each other up to one half the linear dimension ofthe lenslet, 240 microns as in FIG. 8). No doubt there are more thingsto say about this array, but we shall stick with these for now.

Steps 522A&B and step 524 are steps which attempt to address theseissues all in the service of attaining better overall system performance(why else?). Step 522 A&B are philosophically the same and are optionalboth separately and together. They sandwich step 524 since variousnormal considerations enter into whether they should come before orafter or both. It is anticipated, however, that one or the other atleast will be used in normal systems utilizing the invention. As withany set of spatial data, the signal to noise properties of the varioustwo dimensional spatial frequencies of complex₋₋ [][] can be determined,and a "Wiener-like deconvolution" performed on said set of spatial data.Furthermore, it is a matter of normal data processing that such spatialdata can be placed onto another spatial grid ("reinterpolated" as step524 puts it) with minimal adverse consequences, provided the basicsignal to noise considerations are heeded. This is a very general way ofsaying that steps 522A&B plus step 524 resample "complex₋₋ phase []["onto a new "standardized grid" or basis functions, here named"standard₋₋ phase [][][][]" 526, while at the same time minimizing theerror between the "true" complex phase and the estimate on thestandardized grid (i.e. the "Wiener-like" solution) standard₋₋ phase[. .. ] 526 is complex. The reason for using a standardized grid for eachand every field angle derives from later stages where the complex pupilfunction itself becomes conceived as a uniformly changing (continuous)entity across field angles, necessitating a uniform treatment of thebasis functions which describe that function. The reason for minimizingthe estimate error is self-evident. Finally, the "next" statement 502Bof loop 502 is encountered in FIG. 19, with the resulting output of theentire loop being an 8 by 8 array of these standardized complex pupilfunction estimates. In effect, loop 502 merely isolates each field anglein the array of 8 by 8 and performs, generally speaking, prior artmethods on those field angles to obtain the usual phase estimate. Assuch, it is fair to say that all of FIG. 19 is merely applying commonsense adjustments to prior art techniques.

Loop 502 effectively accomplishes the first of the three tasks for theglobal wavefront sub-routine as delineated above. The next task is togenerate (estimate, to use precise language) the corresponding pointspread function (or is UV plane conjugate) which result at the centerfield angle of each and every data sub-zone 320, FIG. 13. It should bepointed out that there are other options at this point, this one merelyis an excellent engineering choice since it is a straightforward problemto convert complex pupil functions into point spread function estimates.It is also compatible with all of the preferred embodiment choices forthe methodology of the local wavefront sub-routine. This issue will befurther discussed in the paragraphs which describe the local wavefrontsub-routines and their advanced treatments. So in order to accomplishtask number 2 of the global wavefront sub-routine, loop 528 is enteredat step 528A, FIG. 19, stepping through all 32 by 32 data sub-zones 320,FIG. 13. Step 530 then interpolates the value of the complex pupilfunction at the center of the current data sub-zone, using the 8 by 8grid of standardized pupil functions 526 as table points. The manner ofthis interpolation is left as a normal exercise in data processing,where it should be noted once again that conceiving of this situation asan "interpolation" rather than a "evaluation of a function at aparticular point" does service to ease of explanation and ease ofprogramming. [It is yet another manifestation of the zonal versus modeldebate which those practiced in the art know so well.] FIG. 12 can onceagain assist in an intuitive supplement to the working of loop 528.There, an array of field angles is represented as a quasi-regular grid306 projected on top of the primary detector boundary. The"quasi-regular" means that geometric distortion can be accounted for ifdesired. Each of the grid points 306 on this array has a complex pupilfunction estimate as given by loop 502. The center of the data sub-zone(see FIG. 13 now, and the 322 projection) then has its own complex pupilfunction interpolated from the actual 8 by 8 grid of the example system.Back to FIG. 19, then, the output of step 530 is a complex array calledlocal₋₋ pupil 532, where once again, this array is standardized and hasan implicit aperture function support. The reason for standardizing thepupil function in step 524 should be clear after considering step 530and considering that other routines can be better optimized by knowing aspecific format is being employed. Step 534 is phrased in FIG. 19 as a"C language" type callable sub-routine named "find₋₋ psf(local₋₋ pupil)"534 where the passing conventions of the C language are also utilized.This routine takes the local pupil function estimate 532, and also using"global variables" of wavelength used (500 nm) and aperture support(circular/annular/etc.), creates a point spread function estimate forthe current data sub-zone. FIG. 19B has a brief and highly generalflowchart for this sub-routine, where those practiced in that art willrecognize this as standard fare: i.e., resampling the pupil functiononto a 128 complex grid, an auto-correlation onto a 256 by 256 complexgrid, fourier transformed into the complex spatial domain, thereconverted to an intensity pattern by calculating the complex magnitudeof the spatial domain complex distribution. Alternatives to this routineare brute force ray tracing, or as alluded to in FIG. 19B, sampling tohigher grid spacings then downsampling back to the 256 by 256 outputarray psf. The output of the sub-routine 534 is a non-negative all-realvalued array named "local₋₋ psf[][][][]" 536, being the estimate of thepoint spread function, at 500 nm, for the field angle corresponding tothe center of the current data zone defined by loop 528. The first twoarray brackets can be seen as indices on the data sub-zone, the secondtwo as the 256 by 256 region. It is a good idea to perform theseprevious operations in floating point precision. The "next" step 528B ofloop 528 is then encountered, so that by the time the entire loop 528 iscompleted, all field angles corresponding to the centers of the datasub-zones now have a corresponding point spread function estimateapplicable to the primary CCD 8 data. This completes the secondappointed task for the global wavefront sub-routine.

The third appointed task, that of further processing the array of pointspread functions 536 into forms which the inner loop 366, FIG. 14, canbest make use of, is thoroughly contingent on what processingmethodology is chosen for creating the specific rev₋₋ array 378, FIG.14, chosen in step 227, FIG. 10. As such, this third task will bediscussed in the context of the third procedure in FIG. 17, which willnow be described.

The best place to begin describing the preferred embodiments of thegeneration of the specific rev₋₋ array 378 (FIG. 14) by the localwavefront sub-routine 376 is depicted in FIG. 20. There shown is ahighly simplified graph, with legend, of a common situation found in theengineering design and in the operation of instruments: the existence ofa multitude of solution to a given problem which variously trade offperformance against resources (or cost). This disclosure presents fourdistinct solutions among a wider array of hypothetical solutions to thespecific reversion methodology utilized inside the main program 350,FIG. 14. [The solutions marked 1A through 1C will be discussed in theparagraphs discussing solution 1]. The graph is deliberately without ascale or without detailed description of the variables on the axes. Itis an intuitive diagram. The four procedures will be in turn discussed.Both solution 1 and solution 2 will be described in detail. The outlinesof the third solution will be presented, not a recipe, and thisdisclosure anticipates that standard engineering development of thisthird solution may turn it into the preferred form of the invention,though at the application data of this disclosure it only exists as anoutline for high priority R&D. Finally, the fourth solution essentiallyentails the "brute force" solution of equation 11, assisted by thegeneral banding (i.e. sparseness) of the H matrix and using thesimplification of equation 13 to best advantage. This solution is welland fine, except that it highly increases the memory, storage, andprocessing time involved. Again, it is anticipated that solution 3 willeventually become "arbitrarily close" to solution 4, primarily due toits formation in the full matrix context (as will be discussed).Summarizing, then, solution 1 will be seen as a very fast method but theone most prone to error, solution 2 is currently preferred, solution 3appears to have the most promise in the near future, and solution 4 isthe brute force full solution that supercomputer users are at leisure toutilize, though at ever diminishing returns.

As an important but ancillary mathematical digression, solutions 1through 4 can be simply related to equations 10 through 13 in thesection on the mathematical background. Future developers of refinementsto the present invention are highly encouraged to take the time toexamine these relationships, since they are intimately connected tooptimizing solution 3. Essentially, solution 1 ignores both thesummation and the convolution in equation 10, effectively treating theproblem as an isoplanatic one, that is, retaining only the diagonalelement of matrix H and setting those diagonal elements equal to the(local) transfer function as defined by psf[][][][] 536, FIG. 19.Solution 2 also ignores the convolution in equation 10, but it does usethe summaries as indexed by q to derive the diagonal elements of thelocal transfer function. Thus, this is still "isoplanatic," but it usesa closer estimate to the effective local transfer function, i.e. therows of H rather than the columns. Looking at solution 1 and 2 in thisfashion can better explain why solution 3 is expected to be preferred inthe evolving future. Solution 3 will merely begin to utilize ever more"neighboring" values of both the summation AND the (heretofore ignored)convolution operation in equation 10. The most elegant intuitive way tovisualize this (see FIG. 25, actually) approach is to first realize thatequation 10 is itself a "shift varying convolution"operation; in thecontext of determining an "effective" inverse transformation for thedata subzone region (i.e., by using a unit g vector, essentially), asimple "deconvolution" can be applied as opposed to no deconvolution atall as in solution 1 and 2, the result still being a better estimateddiagonal matrix. This topic will be picked up in the paragraphs whichdiscuss solution 3. Finally, solution 4 merely solves equation 11 withthe unit g vector (all zeros with a single one) corresponding to thecenter of the data sub-zone.]

I will first describe the simplest and fastest reversions methodologymarked as solution 1 in FIG. 20. This solution is named the "pointspread function" method since it utilizes the array psf[][][][] 536 ofFIG. 14 essentially as is, except that is uses its conjugate. Theprocessing steps are listed in FIG. 21, where those familiar with theart will see that the local wavefront subroutine 376 FIG. 14 is nothingmore than a pointer assignment operation to a pre-calculated (e.g.stored on disk) array of "rev₋₋ arrays." FIG. 21 shows that, for thisPSF-based reversion methodology, the sub-routine 538 in FIG. 19 firstallocates and clears a complex array "global₋₋ rev₋₋ array 549 whichwill be the ultimate output of the sub-routine 538 (be it on disk,memory, or whatever). It then begins with loop 550 through the datasub-zones. The first step 552 inside this loop performs a 2D fft on thearray psf corresponding to the current data zone, where the output goesinto the current global₋₋ rev₋₋ array (real and imag). Thereupon anested loop 554 is encountered which steps through all UV plane points.The processing step 556 inside the nested loop 554 effectively assignsthe new array global₋₋ array[][][][] with the complex conjugatesglobal₋₋ rev₋₋ array 549 by applying the negative value to the imaginarycomponents. Then the two next statements 554B and 550B are encounteredrespectively. In plain English, the sub-routine 538 loads a global arrayof complex conjugated UV plane estimates of the local point spreadfunctions. Finally, the local wavefront sub-routine 376, FIG. 14,reduces to a simple pointer assignment operation (as is quite well knownin software engineering) into the pre-calculated global arrayconstructed to sub-routine 538. For the sake of completeness, this isindicated in FIG. 21A, labelled as step 558.

At this point it is worthwhile to compare this most simple reversionsolution (number 1 of FIG. 20, just described) with prior art techniqueswhich function under isoplanatic assumptions. Reference is again made inparticular to the Primot et. al. paper, previously cited, as an outlinefor the state of the art in this regard. Now it has been "suggested" tome by numerous (independent) people that existing isoplanatic methodscould just as well be "stepped and repeated" across all field anglesmeasured (there still is the issue of how these field angles aremeasured, but in any event . . . ). Since this suggestion has been madeeven by what could be called "non-experts" in the art, I assume thatsuch an idea could be considered as "an obvious extension of the priorart." There are two answers that I generally have for thesesuggestions: 1) yes, it certainly can be used, but it is quite a poorperformer in the error vs. resources graph of FIG. 20, there depicted assolution 1A; and 2) It is quite non-obvious how to move this performancefrom the position as 1A to the position at 1 in graph 20, and I am awareof no efforts, published or otherwise, which have attempted to move froma simple, obvious "step and repeat" methodology to a more sophisticatedone. In fact, I am aware of no attempts to systematically compensate forimagery arbitrarily beyond the isoplanatic region. For the sake ofcompleteness, this disclosure will outline the considerations for movingfrom the "obvious" performance at 1A to the involved reasonings whichhave led to the creation of solution number 1 described above. This"moving from 1A to 1 will itself involve the novel processing methods 1Band 1C referred top in FIG. 20.

FIG. 22 contains a schematic series of methodologies which, moving fromthe top drawing to the bottom, transition from the "obvious step andrepeat" method of isoplanatic prior art processing to the preferredembodiment solution 1, FIG. 20, of this disclosure. We can imagine thatthe performance points of the two intermediate methodologies 1B and 1Ctravelling downward in error but outward in resources as in FIG. 20. Itshould be explained that the reason 1A, 1B and 1C are all out further onthe resources axis is because, in general, they will have a higher noiselevel per acquired exposure, and thus require more exposures to attain agiven level Of noise, or "error" performance. The top drawing of FIG. 22alludes to simply splitting the data region in to an 8 by 8 array ofdata zones each with 128 by 128 pixels each--we are still assuming a 2.5meter telescope, 500 nm, etc. Presumably a wide field wavefront sensoris being employed as described in this disclosure, and thus each regionhas its own estimate of the local point spread function applying to it.The techniques as described in Primo et. al's paper are applied verbatimto each region on a large ensemble of images, and the results are sewnback together into a single 1024 by 1024 final image estimate. It wouldbe fair to say that a post filtering of this final image estimate, wherea spatial filter is applied which attenuates the "square harmonics" ofthe 128 by 128 grid regions, would be part of "the obvious step andrepeat extension of prior art techniques." To reiterate, this "obvious"extension will indeed work and should in general produce results whichare "better than" what would obtain using purely the long term exposureseeing profile. However, whatever average resolution gains are made comeat the expense of gross spatial non-uniformities in signal to noiseproperties (related to the 128 by 128 spacings). Resolution performancewill be non-uniform spatially, and noise properties will be non-uniformin both the spatial domain and the Fourier domain. For visualapplications, central regions can be "degraded" in order to match borderregions, thus at least creating visual uniformity at the expense of aminimal overall gain in resolution. For numeric image science, thesenon-uniformities can cause havoc. Thus summarizing, as mentionedearlier, this "obvious" step and repeat extension of prior arttechniques can indeed improve performance potentially, but only by asmall amount and only within a context which manifests other forms ofperformance degradations.

Stepping out now from extension of the prior art techniques into thenovel design of processing methodologies which borrow on these prior arttechniques, the second schematic in FIG. 22 takes a jump downward inerror performance from number 1A toward number 1's level in the graph ofFIG. 22 (a reminder in again made that this is a highly qualitativetreatment). This technique is there marked as 1B. The idea is to onceagain envision 128 by 128 data regions, perform the techniques of Primotand Rousset to an ensemble of images, only this time only a 32 by 32central "sweet spot" is chosen to be added to a final image. Meanwhile,the 128 by 128 region is stepped over by 32 pixels to the left, the sameprocess performed, and so on, each overlapping region determining itsown 32 by 32 sweet spot. Eventually an entire final image estimate isbuilt on with an outer border of data "thrown out" because it never wasin a sweet spot of any 128 by 128 data region. And once again, postfiltering can attenuate the square harmonics, this time based mainly onthe 32 by 32 sweet spot sewn regions, but also deriving from the larger128 by 128 grid. Looking at FIG. 22 then, experimentation will find thatthe error has subsided at the expense of processing time (however, lesserror can mean less exposures needed, and thus less processing time).Though this new technique 1B clearly gives less overall error, it stillquite computationally redundant and still more error prone relative tothe preferred embodiment of methodology 1, FIG. 20. The second schematicin FIG. 22 also can serve to illustrate a small variation on the method1B just presented. Here, the "sweet spots" are rounded by smoothweighting functions and they are generally broadened beyond the 32 by 32region somewhat (such as bilinear functions extending to a 64 by 64region), greatly minimizing the square harmonics due to the squaresubregions in method 1B, also allowing for smoother "interpolated"values between the centers of the sweet spots. This is a variation whichwas broadly outlined in the disclosure of the preferred embodiment ofthis invention already. Looking once again at FIG. 20, we see a path ofimproving error performance starting with the "obvious" extension ofprior art techniques, moving through novel "borrowings" of prior arttechniques, and then winding up at the superior performance of thepsf-based method of the preferred embodiment of this invention. At thispoint it can be better appreciated the role that the overdetermined setof universal basis functions play in setting up an optimal solution fromfirst principles rather than winding an intuitive path from the obvioussolutions to the optimized. To complete FIG. 22 therefore, themethodology of solution 1, FIG. 20 is sketched as a next logical"intuitive" step where the redundancies of methodology 1C areessentially eliminated (as well as the concomitant error sources) and a(maximally?) efficient psf-based methodology is arrived at. The wholereason for outlining methods 1A through 1C was to show that "obvious"wasn't necessarily so obvious, and when it is indeed obvious, it isgenerally of limited value. That is to say, even extending basicprinciples of prior art techniques requires a fair degree of novelty.

We can turn now to a description of the preferred embodiment of theprocessing methodology; specifically, the use of the data point windowfunction based reversion method, labelled as solution 2 in FIG. 20. Nowthe definition of a"data point window function" was already given in themathematical background. To revisit this definition in plain English, itis the weighted sensitivity function that a given pixel has on objectspace at any given moment in time. More prosaically, it is the rows ofthe H matrix in equation 1 rather than the columns. [My apologies forusing the term "window" since this may be easily confused with its usagein the spectral estimation context; i.e., the weighting of spectralcomponents in order to affect preferable space domain properties.However, I must say that I feel the use of this "mother's tongue" wordis infinitely more appropriate to the manner in which it is presented inthis disclosure than in its more stolid spectral context. It is quitethe matter for aesthetics as well as agreed jargon, but admit, a pixeldoes truly have a "window" on object space, does it not?] The very factthat it is a row rather than a column should give a first clue as to whyit is a better performer than the PSF--it directly relates to the datapoint under question (whose reversion profile is being sought) ratherthan indirectly as in the PSF. As mentioned earlier however, it stilltreats the overall reversion problem as a quasi-isoplanatic one (it canbe likened to having the row define an approximate Toeplitz matrixrather than t he column defining that Toeplitz matrix as in the PSFformation). It will take solution 3, FIG. 20, to take the approximationyet further into the true anisoplanatic realm. But for now, the datapoint window function methodology works admirably, if not perfectlyoptimized. [It is important to none that the very idea behind the datapoint window function posits a space variant matrix H. In the invariantcase, the point spread function and the data point window function aretrivial mirrors of each other and there is no reason to distinguish thetwo. I have always assumed that this is why there has been littleattention paid to t he distinction in the first place. As can be seenhere, however, the distinction is quite important.] The price paid forusing the data point window function is generally one of increasedcomputation time (and usually more complex software data managementschemes). This derives primarily from the fact that the PSF is a simpleconjugate to the complex pupil function along a given field angle,whereas the data point window function generally must be composited froma neighborhood of point spread functions (in practical terms). This canbe seen more readily when we describe the details of the processingsteps as laid out in FIG. 23.

The output of the sub-routine 538, FIG. 19, for the data point windowformulation is generally the same as it was for the PSF formulation,except that the global₋₋ rev₋₋ array is now called the global₋₋ window₋₋array 561, having the data point window functions stored as opposed tothe complex conjugates of the point spread functions. FIG. 23 beginswith a description of the sub-routine 538 which should be performed whenusing the data point window reversion methodology. It first notes theexistence of the previously calculated local₋₋ psf array 536. Theprocessing steps begin with step 560 which allocates and clears a 32 by32 by 256 by 256 array, global₋₋ window₋₋ array 561, of unsignedintegers (floats can be use if memory is no problem). This array willhold the calculated window functions as they exist at the center of thedata sub-zones 320, FIG. 13. [For practical computing purposes, we shallcalculate the window function at the center of the pixel at the 16th rowand the 16th column; rather than the fractional "true" center at(16.5,16.5)] Recall that it is a simple matter to calculate the columnsof the H matrix in equation 1 based on the complex pupil function. Infact, the preferred embodiment ONLY utilizes the calculation of thecolumns (the point spread function) as a means for transforming betweenthe pupil plane estimate and the spatial domain estimate. Yet the use ofthe data point window function necessitates calculating (estimating tobe precise) the rows of the H matrix. Clearly we could calculate eachand every column of the H matrix, i.e., each and every point spreadfunction in the image, and then fill all row estimate sin the doing, buthis requires massive computational time and resources and is ratherinformational redundant. Instead, the preferred embodiment performs theoperation which is described in the second, highly simplified intuitivediagram of FIG. 24A and 24B--it still looks complicated even though issimplified, admittedly. What is there depicted is that 1) We only wishto calculate rows every so often, effectively at short enough intervalsto adequately represent the spatial variance of the H matrix, in ourcase being at the center of 32 by 32 pixel regions (i.e. generallyisoplanatic); 2) We fill all the chosen rows by stepping through columnsat a comparable spacing as the row samples, where we will again use 32by 32 spaced regions; 3) At each columns, we calculate its values (itspoint spread function) according to "find₋₋ psf" of FIG. 19B; 4) thenstepping through the rows which we are calculating; 5) effectivelyrotate the current column 90 degrees clockwise about the common elementof the current row and column, i.e. the "current axis"; 6) then multiplythe rotated values by an "interpolating function" centered about thecurrent axis and add the resultant into the accumulating row values.[This figure is most simplified because it ignores the subtleties of a 2dimensional image-type matrix, i.e. it ignores multiple bandingsubtleties. It also only depicts the column being rotated onto threerows, where in general it will be rotated onto a total of 81 rows (9 by9) as in the lower left drawing of FIG. 24A.] The importantphilosophical point to note is that we are using the slowly changingproperties of H matrix (its slow spatial variance) to chunk it up intosmoothly interpolated quasi-Toeplitz mini-regions, then performingsimple transposition operations that Toeplitz matrices allow. The wholepoint of this fanciness is a) I have found no easier way to directlycalculate the data point window functions, and b) it greatly saves oncomputational time and resources, that is, it is more informationaleconomical.

The remainder of processing steps in FIG. 23 translate the concepts inthe two drawings of FIGS. 24A&B into computational routines. The nextstep after allocating the large array 561 is step 562 which allocatesand "loads" a 64 by 64 floating point array "weight" 563 of valuesrepresenting a bilinear interpolating function. [Smoother interpolatingfunctions than the bilinear, or ones of larger extent, most certainlycan be used instead; usually with a larger than 64 by 64 region.] Usingthe C convention of numbering arrays with the first element as 0, eachelement in the array will have the value:

    bilin[i][j]=(1.0-32i|/32)(1.0-|32-j|/32) (16)

[Efficient software engineering may wish to use quad or octal symmetryto merely store and use a smaller segment; but the computationaloverhead in this routine is so small that using the full array can makefor easier pointer arithmetic.]

Step 564 then creates a new virtual (hypothetical, computationallyspeaking) array 565 of 40 by 40 data sub-zones, the 32 by 32 real datasub-zones plus an extension of 4 zones on all borders, accommodating the128 pixel overshoot of the 256 by 256 universal basis functions. Next,loop 566 is encountered which steps through the indices of this 40 by 40virtual array. The first step 586 inside loop 566 is a query as towhether the given virtual zone is inside or outside the real 32 by 32data sub-zone region; if it is outside, that the closest real datasub-zone's psf array 536, FIG. 19 will be used in the ensuing steps,otherwise, if it is inside, than the existing array psf is used. Thus,step 568 performs a very simple extrapolation if necessary; moresophisticated schemes are again quite employable as well. Step 570 is anested loop which will cycle through the "applicable" window functionsites, i.e. those sites which are in general within the 256 by 256 pixel"reach" of the point spread function. This "applicable" region is ingeneral a 9 by 9 array of window functions sites when the currentvirtual zone defined by loop 566 is well within the real data region,and it is truncated from a 9 by 9 array when the virtual zone comeswithin 4 zones of the border and beyond the border. [This is soundingmore complicated than it is; Those practiced in the art of imageprocessing will recognize this as the usual "trivial" that comes withdealing with data borders.] FIG. 24A's lower left drawing has aschematic helping to better explain the "applicable" criteria. Yetanother nested loop 572 is then encountered which steps through all 64by 64 indices of the weight array, effectively stepping through theregion 575 in FIG. 24A. Inside loop 572 there is first another query 574as to whether a given point is inside the psf region 536 or outside(since the region 575 can be partially outside itself). Then step 576 isan accumulating weighted add, where the indices on the accumulatingresultant, global₋₋ window₋₋ array 561, is mirrored about the centerpixel (128, 128) from that of the psf array. This mirroring is indicatedin the two lower right diagrams of FIG. 24A. [Had we used a simplesquare box region of value 1.0 as our interpolating function instead ofa bilinear, step 578 would have been a simple assign statement.] Itshould be noted that several offsets are involved in these loops andcalculations, and it is no small software engineering task to keep trackof them all. The important point to realize is that a mirroring ofindices does take place, and care should also be taken with theconventionals of using the (N/2)'th element or the (N/2-1)'th point forclarity. Three "next" statements are then encountered, closing loops572, 570, and 566 respectively. Step 577 then indicates that ALL of the32 by 32 array of window functions must then be fourier transformed intotheir UV plane counterparts, using the same 2D FFT routines mentioned inthe main program of FIG. 14. The result is then a filled array "global₋₋window₋₋ array" 561 which is merely the global set of (complex, UVplane) rev₋₋ arrays for all the data sub-zones. Hence, sub-routine 376,FIG. 14 becomes a simple assign loop (or pointer assignment in the "C"language), represented as step 580 in FIG. 23, just as in FIG. 21 in thepsf method. Actual hardware systems and software programs may find thatthe array sizes involved are so large that a fair amount of "swapping"is needed, or, some of the steps in the sub-routine 538 of FIG. 23should be done inside loop 366, FIG. 14 rather than outside. This is whythere is the sub-routine called the "local wavefront sub-routine 376,FIG. 14," even though for solution 1 and for solution 2 of FIG. 20, bothhave turned out to be simple pointer assignments.

This completes the detailed description of the preferred embodiment, theone which uses the data point window function as the basis for thereversion methodology. 1 turn now to a discussion of what is the broaderde facto preferred embodiment of the invention. It is "de facto", sinceboth of the preceding solutions 1 and 2, FIG. 20, are effectivelycontained within its scope.

Recalling the discussion surrounding equations 10 through 13, andturning to the drawings in FIG. 25, we can outline the technique forrefining the specific reversion methodology to be used. The top drawingin FIG. 25 presents a simplified view of a data point centered at theorigin for convenience, having the point spread function of the nominalfield angle of that data point drawn an x with a circle around itdirectly at the origin as well. Likewise, there are field angles nearbywhich have slightly different point spread functions, where the distancescale here might be spaced 1" apart, somewhere within its isoplanaticregion which would isolate the "low order" variability in the pointspread function. Those practiced in the art will realize that there isstill a high degree of correlation between the point spread functions atthese nine locations, but in some sense there is a local "tip-tilt andlow order curvature" present in the field variant point spread function.In other words, if we focus in on only the local low order change in thepoint spread function in the neighborhood of the field angle at theorigin. We can therefore imagine that modelling "low order" change inthe point spread function in this small neighborhood is better thanmodelling no change at all (as in the case of solution 1 and 2, FIG.20). The goal now becomes the incorporation of a finer derivation ofrev₋₋ array 378, FIG. 14, based on that local change. It appears thatthe best way to approach this problem is to view the localcharacteristics of this data point region. To do this, the lower drawingcontains a highly simplified matrix equation which attempts to"localize" the situation depicted in the top drawing. The text underthis figure further explains the approach toward a solution 3, FIG. 20.It will be found that by assuming that the point spread functions of the8 neighboring field angles (to the central field angle) are extrapolatedinto an arbitrary distance, the matrix equation 11 simplifies into thatdepicted in the lower drawing of FIG. 25. By placing an all 1's vectorin Gn. we set up a solution for the inverse (as opposed to "reverse")projection from the data point at the origin into the solution space.Solution 1 and 2, FIG. 20, effectively did this same operation, exceptthat only the diagonal of the H matrix was used, not a small hand as inthe figure. So without going into any further detail, it is anticipatedthat further refinements will be made which determine yet finersolutions than the data point window function of FIG. 23.

Finally, solution 4 of FIG. 20 is merely cranking out the full (huge)matrix solutions to equation 1 and/or 11. There is truly no magic inthis, it is an a priori given that it can be done. However, it isgenerally considered too bulky and inefficient in any event, andbesides, it does not pose the multi-data set reconstruction problem inthe overdetermined least squares formalism of FIG. 14, thus "arbitrary"methods of dealing with nulls must be used.

Advanced Wavefront Estimation Routines

Whereas the processing methodology presented in FIG. 19 and discussed inthe disclosure text was quite functional and "enabling," it cannot besaid to be the "best mode" of the invention. The only reason for this isthat it tends to have higher residual wavefront error that the routineswhich follow. Referring back to FIG. 5, there is box 41 which speaks of"pseudo-inverting" the spatial filtering which is taking place opticallyat the focal planes 66 and 98, FIG. 8. Loop 506 and step 518, FIG. 19,in particular is where our focus will once again turn for improvement,though not necessarily inside the loop, but in modifying it after it andthe traditional wavefront estimation routine of step 518 is completed.In loop 506 and step 518, the procedure essentially followed basic priorart Hartmann methods in determining the wavefront estimate and called itgood. In the prior art systems, however, spatial filtering was much lessof a consideration and/or error source. When the wide field designs, theissue becomes central.

FIG. 26A gives a closer look at what is taking place inside the widefield Hartmann camera system. A bright star 601, slightly offset fromthe nominal central field angle of one of the 8 by 8 field Hartmannsampling apertures (e.g., 46, FIG. 4) is producing a polychromaticspeckle distribution at the nominal focal plane 66, FIG. 4 (thisaperture would be a square in the example system of FIG. 8). Thespectral distribution could be something like the spectral curves ofFIG. 12A, right diagram. The sampling aperture 600A of FIG. 26A is anidealized sharp edged profile of a sub-aperture, where inside the circleit is transmissive and outside it is completely absorbing. Thepolychromatic light distribution (still complex, i.e. this is NOT anintensity profile, per se) is drawn overlaid in iso-modulus lines, wherethose area of the light field transmitting are drawn as solid lines andthose absorbed are drawn as dashed lines. Modern optical language andFourier Optics theory therefore tells us that the sub-apertures areperforming a spatial filtering operation on the two conjugated planes:the telescope pupil plane 34, FIG. 2 and the lenslet sub-aperture plane52, FIG. 4. The drawing of the bottom of FIG. 26A shows that thefiltering function in actuality will be a slightly rounded distribution,a fact which is important when considering the 2X and Full wide fieldHartmann designs of FIGS. 6 and 7 respectively, since these designs inparticular have 45 degree inclined apertures giving rise to somewhatasymmetric "filtering functions" 600B. We also note that spatialfiltering is a decided function of wavelength. Those practiced in theart will recognize that the complex pupil plane wavefront has beentransformed into the wavefront existing at the lenslet plane 52 via aconvolution operation defined by the effective aperture function drawnas 600B. The spatial offset of the bright star 601 manifests itself asthe usual phase shift (see bottom equation, FIG. 26B). The result of allthis general theorizing is the realization that the measured values ofthe wavefront slope are inherently altered from the slopes existing atthe true telescope pupil plane 34, FIG. 2. Yet it is the true pupilplane function we are seeking for the processing methodologies, andtherefore we need to seek methods which better estimate the true pupilplane function rather than making do with the measured values as in thenominal procedure of FIG. 19. Those familiar with the art of fourieroptics should immediately recognize that, indeed, better estimates ofthe pupil function can be found, but within other sets of limitationsand constraints; they will also realize that there is a range of toolswhich can be brought to bear on creating better estimates rather thanone singular "method." The remainder of this section therefore discussessuch solutions in detail, including their philosophies. It should alwaysbe borne in mind that the end goal of all of this is minimizing theresidual wavefront error between the wavefront estimate created by theprocessing methodology and the "true" wavefront.

FIG. 26B sets up some generic mathematical concepts and a frameworkwhich can be used in creating computational solutions to improving thewavefront measurement/estimation. The Full or 4X wide field Hartmanndesign of FIG. 7 and FIG. 8 will be the basic model used in the drawing,where the 3.6" field angle spacing of the specific preferred embodimentis also used. The top left drawing in FIG. 16B is a three-dimensionalarray representing the "current" multi-spectral object estimate 610corresponding to one of the field angles in the 8 by 8 array of sampledfield angles. It is a course resolution estimate of 1" square boxes of4" square total extent, giving a 4 by 4 pixel estimate. The doublesummation below merely tracks the coefficients of the values and leavesoff the implied square box basis functions, as is common. The spectralplane can be considered as 20 nanometer square bandpass regionscorresponding generally to the spectral regions of significant energy asdefined by the rightmost spectral curves of FIG. 12A, the wavefront beamsensitivity function. If a monochromatic "artificial beacon" is the onlysignal of significant energy content, then there will only need to beone spectral plane rather than several, where this beacon assistedsystem is also depicted in FIG. 12A--efficiently reducing it to a 2dimensional object estimate. The course 1" square boxes of 4" squaretotal extent, giving a 4 by 4 pixel estimate. The double summation belowmerely tracks the coefficients of the values and leaves off the impliedsquare box basis functions, as is common. The spectral plane can beconsidered as 20 nanometer square bandpass regions correspondinggenerally to the spectral regions of significant energy as defined bythe rightmost spectral curves of FIG. 12A, the wavefront beamsensitivity function. If a monochromatic "artificial beacon" is the onlysignal of significant energy content, then there will only need to beone spectral plane rather than several, where this beacon assistedsystem is also depicted in FIG. 12A--effectively reducing it so a 2dimensional object estimate. The course 1" resolution of "the estimatepixels" derives from several considerations. Perhaps the larger reasonis that until the whole technique of wide field wavefront sciencebecomes a meticulous, mature science, this may be all the resolutionwhich is required to effect "reasonable" error reduction. Besides, theresolution certainly can be increased as will be discussed later. The 1"resolution of the Figure also implies that an ensemble of the primaryspeckle image themselves, filtered by gaussian distributions forinstance, can serve as the initial object estimates. A "white light"model can also be used as an initial estimate, where all spectral planesthen have identical values. These broader issues will be discussedimmediately following the description of FIG. 26B, especially in FIG.28. Turning now to the upper right drawing in FIG. 26B, there isdepicted iso-radiometric lines--10%, 50% and 90%--for the amount ofoptical energy that the in situ aperture receives as a whole, laid outon the same spatial reference of the object plane. The 4" square box ofthe object estimate is drawn as a reference. The effects and subtletiesof the 45 degree inclined apertures depicted in FIGS. 6A and 6B arealluded to here. The fact that one edge of the aperture is farther outof focus than the other (the effect is not huge by any means) is drawnas a smoothing and broadening of the iso-radiometric lines to the rightof the figure. For reflective telescopic systems, the radiometricprofile should be effectively the same across wavelength, and if itisn't, and if the difference produces appreciable error effects lateron, then these radiometric profiles certainly can be treated as amulti-spectral entity as well. On the other hand, when the profile inthe upper right drawing of FIG. 26B is fourier transformed into W andmodelled inside the "convolution" equation at the bottom of FIG. 26B,then there will be a scaling factor which is directly related towavelength, hence the lambda subscript on W. The equation at the bottomtherefore describes a mathematical model of the operations performed onthe "true" pupil function, there named "pupil," and the measured phaseestimate. It can be seen that each "spectral pixel" has its own uniqueshift and convolution applied to it. The pupil function is taken as aninvariant function in this figure and further refinements can begin tolook at its local tip-tilt and low order curvature (in that event, the"pupil" function would thereby be given the subscript in denoting thatit changes as a function of field angle). The s_(m) is a two dimensionalspatial coordinate representing the shift of the pixel from the nominalcentral field angle. The trick now becomes how to turn this mathematicalmodel of the degradation into methods which better estimate the actualpupil function.

Among the many choices of better estimating the true pupil function, thepreferred embodiment of the invention chooses a rather simpleapproximate scheme loosely based on what has generally been known as theVan Cittert solution to deconvolution problems. [The basic methods setdown by Herr Van Cittert in the late twenties and early thirties havegone through the nth degree of refinement and modification.Nevertheless, using the "Van Cittert" label still carries the basicphilosophy of solution with it.] It is a little bit ironic that thisdisclosure chooses such a solution, since much of the philosophy of theinvention as a whole is to boil down certain relationships into direct(quasi-deterministic) forms rather than the more indirect Van Cittertformation. Be this as it may, the Van Cittert solution will indeedeffect improvements at a reasonable computational cost, and futureengineering refinement will be tasked with finding betterquasi-deterministic methodologies which are better candidates--as ifthey were plotted on the same type graph of FIG. 20.

FIG. 27 has a general processing flowchart of a new sub-routine whichshould be placed immediately following step 518 in FIG. 19. It istherefore called sub-routine "wavefront₋₋ est₋₋ refine" 519 as if itcomes right after that step. The processing methodology begins with theusual definition of what it expects as input, in this case being thecritical entity of the phase estimate at the pupil plane, as derivedfrom traditional wavefront estimation routines in step 518, FIG. 19. Thefirst processing step 620 is allocating and clearing two new arraysnew₋₋ phase[][] 622 and tmp₋₋ phase[][] 623, each with the same numberof elements as the array "phase" which is the output of step 518, FIG.19. The next step 624 is a step which need only be performed one permain program, where two "table arrays." w₋₋ lambda 625 and e₋₋ shift626, are pre-loaded with vales corresponding to the fourier transform ofthe aperture function and values of the phase shift of a given(spectral)pixel, corresponding to the situation as described in FIG.16B. If necessary, the w₋₋ lambda array could be a further array of the8 by 8 field angles if each of their aperture functions is sufficientlydifferent. Step 628 is a sub-routine which loads the values of the lowresolution spectral object estimate of the top left drawing in FIG. 16B.Recall that this can be a gaussian filtered (FWHM 1") of the ensemble ofprimary speckle images, then each field angle assigned its own unique 4by 4 pixel region; the spectral characteristics can be treated asuniform if no explicit information is given to the contrary, and if itis treated as "spectrally neutral," only one array f_(m) need becalculated. Then we encounter loops 630 and 632, stepping through lambdaand then through the low resolution pixels respectively. The order couldbe changed here if thee are software efficiency issues at stake. Thenext step 634 is complex operation multiplication on the named arrays.Here, the array tmp₋₋ phase is assigned a phase shifted version of theinput phase array which was the original output of step 518, FIG. 19.The phase shift array 626 was pre-calculated and is specific to eachlambda and each pixel indexed by m. Step 626 is also a complex operationstatement, where now the array new₋₋ phase created and cleared in step630 has an accumulating add (accumulating inside the lambda.m loops) ofthe product of the scalar value f_(ml) multiplying the convolution ofthe tmp₋₋ phase 623 loaded in step 634 and the table array w₋₋ lambda625. Then the "next" statements of loops 632 and 630 are encounteredrespectively. Then a new loop 638 is begun, where the indices on thephase arrays are explicitly steeped through rather than implicity. Theinner statement is then the crux of a Van Cittert like operation, wherethe difference between the incoming phase and the newly calculated phaseis found, and some empirical multiplication factor alpha is applied tothe difference, and the difference duly subtracted from the originalinput array "phase." Using an alpha of somewhere between 0.5 and 1.0 isa crude rule of thumb; empirical values always seem to wind up being thefinal choice however. Finally, loop 638 is closed with the "next"statement 638B. The output of the sub-routine is the same as the input,in other words, only with than the raw slope measurements would haveotherwise given.

Again, it should be clear to those practiced in the art of imageprocessing and in general "data correction methodologies" that theexample processing methodology of FIG. 27 is indeed indicative of anentire genre, iterative and quasi-deterministic alike. Entire books havebeen written on the basic subject and these building on the present workare encouraged to consult those books and apply the basic solutions tothe specific conditions at hand. With this said, we can turn to a briefdiscussion of FIG. 28 wherein a rough refinement feedback system isdepicted on the processing methodology of the invention as a whole. Atits core, this could be conceived of as a grand web of Van Cittert likefeedback loops which in principle can iterate its way to a refinedself-consistent solution for an object under study. However, I believethis interpretation overstates the argument for the very need for theserefinements, and in situ testing of actual telescopic systems will bethe final arbiter of whether such loops are employed. Furthermore, themore that such "feedback loops" are replaced by direct, efficientprediction methods, the better. In some senses, the Van Cittert approachis hence a fail safe when there is appreciable error which can beremoved, but the quasi-deterministic solutions are eithercomputationally impractical or they are not within the grasp of theimagination. In any event, some or all of the loops depicted in FIG. 28should be considered as a potential part of a working model of theinvention. I personally like to envision them as single iteration loops,since the errors they are correcting are already very small (except inthe case of FIG. 27 perhaps), and one iteration will bring them wellwithin the point of diminishing returns. The elements of FIG. 28 thatshould be particularly noted here are first the step labelled "spotcross-talk correction" 704. Very briefly, it can be imagined that afterall spots of all field angles have been found and centroided, a certaincross-talk predictability becomes apparent which will tend to shift thecentroids of two neighboring spots slightly closer to each other as thetwo spots become closer themselves. That is, the positive slopes of aneighboring spot skew the 9 point centroid in that direction (see FIG.18A). As one spot moves in the direction of one neighbor and away fromothers, the effect becomes pronounced and asymmetric. Provided with acrude model of spot light distribution, a small correction factor couldbe applied to the centroid values of all the spots. The real importanceof refining such a scheme as this is that, ultimately, a finer fieldangle spacing of sampled complex pupil functions could be managed (say3" spacing instead of 3.6") if and when the spots can be more closelypacked together. Spot wander and outright overlap provides a harderlimit on how closely packed the field angle sampling (and the spots) canbe. Step 706 poses the possibility of performing an operation such asdescribed in FIG. 27 on the Hartmann digital images directly, ratherthan on their derived phase estimates. More broadly, local low ordercurvature of the phase function can also be incorporated. The resultingrefinement procedure would use the error between actual spot data and"predicted" spot data as the modification signal on t he pupil functionestimate itself. This may be time consuming, or it may not be, dependingon how much computing resources can be brought to bear on the problem.Step 708 is precisely what was described in FIG. 27. It had its ownFigure and its own highlighted explanation since it is anticipated thata large majority of applications and telescopic system error budgetscould use its error reduction magnitude. This general question goes forall such "refinement schemes" obviously. Step 712 asks whether after anentire solution 412, FIG. 14 is found, would a better solution be foundif this estimate is used as the object estimator in the top left handcorner of FIG. 26B in the wavefront sensing routine. Step 512 isspeculative (in terms of whether it is past the point of diminishingreturns) at this point, but its possibility certainly exists. Finally,another speculative option is considered in step 714, suggesting that inapplications where a spectral set of images is being taken of a givenobject, would this information be sufficiently helpful in improving thespectral plane estimates of the upper left hand drawing in FIG. 26B aswell. Again, whether or not this could substantially improve overallperformance is not the issue here, only the existence of the possibilityif needed.

A few additional refinement possibilities which have or haven't alreadybeen mentioned and which did not easily fit into FIG. 28 are thoseconnected with improving the generation of "rev₋₋ array" to the ideal ofthe full matrix solution of equation 14 (i.e., yet more discussion onthe search for solution 3, FIG. 20), as well as refinements whichexplore further into the role of "local low order shift variance" andhow such simplified viewpoints can simplify both the reversionmethodologies and the wavefront estimation methodologies. Finally, thescintillation function of the pupil plane as usual is given short shriftand all too often given the a priori boot. Perhaps accounting for it andusing its "measurement schemes" will be well into diminishing returns,perhaps not. The disclosure already had a short essay on this earlier.In any event, rementioning it here in the context of FIG. 28 and the"other" refinements possible is not inappropriate.

Concerning polarization effects and their potential introduction offurther error sources, this is an issue present in this invention and inthe prior art systems. Once the basic solution has been formulated inthe manner of equation 1, the isolated issue of polarization becomes oneof predicting the rows and/or columns of H within this additional"polarization" context. The issue of how the "spot information" on theHartmann data is different from the "true pupil function" and the focalplane array 8 light distribution is of equal concern in the engineeringdesign of an embodiment of this invention as it is in adaptive opticssystems, and should be given the same engineering attention accordingly.

One final note is on the broad issue of tracking an object as the earthrotates. Though there are volumes which can be written about this andits effects on the methods of this invention, perhaps the easiestsummary of the effects of "tracking errors" is that they should, infact, blur the final image in the same manner that they have alwaysdone. At first glance, this appears to place major new demands on thetracking abilities of large telescopes, suggesting they go from therecurrent 0.1" region specifications to an order of magnitude better (Ican hear opto-mechanical engineers groaning). Not to worry however,there are several interesting tools which can be brought to bear onalleviating this problem. POSSIBLY TO THE EXTENT WHERE TRACKINGSPECIFICATIONS CAN ACTUALLY RELAX. Wouldn't that be a twist. The basicidea is that phase only corrected version of a speckle image should havestable overall "centroid" properties, and that by either taking theindividual exposure numerators RE, IM, 384, FIG. 14, or taking these andperforming modest modulus amplification, these entities will exhibitcentroid figures which should change from one exposure to the nextaccording to A) the tracking, and B) random noise from instrumentationand shot noise. The magnitude of the noise effects on the centroidshould be quite tolerable for brighter object (scenes), while for lowlight level objects the situation becomes less useful. So at least tosome "brightness-level," this centroid correction can b e used tocontinually shift the image estimates according to a global tip-tiltparameter applied to the data. For the lower light situations, ratherthan attempting to design absolute rigidities of the components of anoverall telescopic system plus a rigid bright star tracker data feedback(all of which can be used as data in set 218, FIG. 10 to be sure),another option is to merely monitor the orientations of the primary CCD8, the primary mirror, the secondary mirror, and a star tracker veryprecisely, let them be relatively "non-rigid" because of the shortexposures, and be content with merely knowing their orientations ratherthan "enforcing" rigid orientations. This could be done withinterferometric techniques among others. In any event, chasing the lowerbrightness tracking error suppression techniques is yet one more areafor detail engineering.

This concludes the disclosure of the preferred embodiment of theinvention. The following sections describe major variations orextensions thereof applied to systems which are different enough todeserve being highlighted.

Alternative: No wide field wavefront distortion information available

Modified methods of the invention can be applied to physical systems andto sets of blurred digital images wherein no direct information isavailable concerning the distortion of the images. In effect, thesituation would be such that there would be no input information of thetype 216, FIG. 7. Such a situation can arise for many reasons, and thepast two decades have generally referred to this problem as "blinddeconvolution" among other things. In the case of this invention, andnot withstanding the upcoming description of how the techniques of thisinvention can be employed to improve the performance of conventional"isoplanatic systems" the teachings of the present disclosure can beapplied to the restoration of blurred digital images which altogetherconstitute a more general class of "field variant blind deconvolution."The distinction is important since there are published prior arttechniques which explore rather deeply the issues and solutionsapplicable to "field invariant blind deconvolution" (ref. Bates, &McDonnell, "Image restoration and reconstruction," The OxfordEngineering Science Series, 16, Clarendon Press, Oxford, Section 10.),i.e. isoplanatic solutions. Even so, the methods of this invention maywell present superior results to imaging deblurring applications whichare forced to make isoplanatic assumptions (field invariance) even inanisoplanatic systems.

Getting back to the issue at hand, the case we examine now is when thereis no input data set 216, FIG. 10, available. In general, it can beassumed that significant system data information as well as grossstatistical information about the distortion functions IS available asinput 218, FIG. 10. And of course, the set of "E" distorted digitalimages must be part of the input data 214, FIG. 10. For the sake ofsimplicity, we will again assume that a 2.5 meter telescope is providingthe distorted images and that an 8 by 8 field of "unknown distortionestimates" would be sufficient to characterize the distortions. Eachdistorted digital image in the input set 214 sets up its own"Z-dimension" probability region in solution space, where thatprobability is determined largely by the gross statistics of thedistortions contained in input data set 218. Since Z is effectively<<N², the number of data points deriving from the primary detector, aunique high probability region should emerge in solution space afterseveral distorted images have their respective probability functionsmultiplied together. It is not the point of this disclosure to furtherexplore the subtleties of t he higher mathematical reasonings involvedwith this problem. Instead, it is rather clear that a simple methodologyfor searching for and finding this region of higher probability, i.e.,the sought solution 224, FIG. 10, could be facilitated by the techniquesof this invention. The two practical difficulties which will attend anysystem incorporating these methodologies will be, in general, the muchlarger computational times involved and the traditional nasties involvedwith search algorithms (slow convergence and local maxima to name thetwo largest). For solutions to and/or help for these common searchalgorithms, there are many texts and published accounts available. Also,there are certain basic limitations to the results of these techniqueswhich are not only brought about by basic signal to noise ratioconsiderations, but also due to certain quasi-non-uniqueness which (inatmospheric turbulence at least) exist within the gross statistics ifthe distortion functions themselves. An intuitive situation is where thezernike coefficient connected with the focus adjustment, though it isonly one dimensional, can blur a given image only along a unique path insolution space; several manifestations of "defocus" (alone) thereforegive rise two mostly similar families of probability distributions,leaving little assistance in searching for a unique "intersection point"or high probability point. There are pragmatic solutions for items suchas this, but they do tend to limit the results obtainable. In any event,though this degenerate case of having no direct information 216, FIG.10, is far from the preferred embodiment of the invention, some systemsmay have no other choice, and the processing techniques described fromFIG. 10 onward may not completely "solve" these problems in light of theproblematic search difficulties just mentioned, but they will no doubtprovide significant tools in such bona fide "invention."

Special Section: The use of artificial beacons

It has long been appreciated that wavefront sensors generally require amoderate incoming light flux in order to operate properly. For the caseof atmospheric turbulence, their short exposure times increase the needfor higher light levels. Both theoretical limits and practical state ofthe art instrumentation limits have generally limited conventionalwavefront sensors to the 10th stellar magnitude object brightness orperhaps a little better. [It is actually dependent on a variety offactors.] Within the last decade or so, practical systems have beenbuilt which provide additional light levels to the wavefront sensorsvia, generally, laser light directed along a near-paraxial field angleof the telescope; the laser light either backscatters back into thetelescope, or may instead excite fluorescence in higher altitudeatmospheric layers with a small portion of the energy once again goingback to the telescope. These systems now go under the name of"artificial beacons" generally. There ultimate purpose in to extend thefaintness limit of objects which can be imaged by the telescope (almostexclusively in the context of adaptive optics systems; see thediscussion on prior art). These artificial beacon concepts can also bequite helpful in conjunction wit the techniques of this invention.

The wide field nature of the preferred embodiment of this inventiondictates that artificial beacons should cover the field of view desired,as opposed to "some close location to an object of interest within theisoplanatic region", as is often currently the case. This suggests a64-fold increase in optical poser (not a trivial matter) from the sourceof the artificial beacon over and above that required for theisoplanatic case, merely deriving from the 8 by 8 isoplanatic regionsbeing covered (and assuming comparable limiting magnitudes of theHartmann sensors used). We have already seen how the wide field Hartmannsensor integrates its signal from an extended field of view, thus, thereis no inherent reason to localize the beacon light into a "star." Ifparticular applications could benefit from each localization, so be it.The point is that for wide field uses of artificial beacons, a very wideartificial beam can be generated and projected out along the field ofview of the detector: there is no need to get any fancier. Suffice it tosay that the techniques and methods of this invention can certainly makeuse of this latest methodology of "artificial beacons."

Special Section: Non-Imaging Uses of the Wide Field Wavefront SensingSystem

It is well known in the optical industry that wavefront sensing systemscan be employed as central or peripheral elements in non-imagingapplications. It is clear that the wide field wavefront sensing systemdisclosed herein is no exception. Specifically, the disclosed systemcould certainly be utilized in applications such as single or multi-optical surface metrology (profiling, etc.), distorting mediatomography, electromagnetic communications systems including laser-basedcommunications through a distorting media such as the atmosphere, and soon. Overall, it is quite likely that the informational basis ofanisoplanatic wavefronts, as measured by and provided by the disclosedsystems, will offer significant performance increases, or outrightenablement where no enablement was possible before, of many of thesenon-imaging uses of wavefront sensing systems.

The Wide Field Hartmann Shack Wavefront Camera

This section describes an improvement to the full field wavefront sensordescribed above. It shall simply be referred to as the Wide FieldHartmann Shack (Wavefront) Camera. It has a less complicated hardwareand optical layout, it is simpler on a conceptual plane leading to asimpler approach to the processing methodology, and it should be able tofunction on dimmer objects. As was the case with the above-describedfull field wavefront sensor, this wide field Hartmann Shack (HS) camerawill function on extended objects as well as on arbitrary contrastextended objects (including zero-constrast "blank walls"). The onlypotential drawback of the wide field H/S camera is that it may haveslightly less immunity to non-singular solutions in every highturbulence conditions relative to the beam-splitting full field design;that is, there may be less tools available to "disentangle" overlappedand inverted signals on the wavefront distortion. This latter potentialdrawback may be quite minor and may only be of consequence in extremelyhigh turbulence applications not typically encountered in astronomicalsituations.

The wide field Hartmann-Shack camera is an extension to the generalphilosophy of the traditional Hartmann-Shack wavefront sensor. Whereasthe traditional H/S sensor attempts to measure the local tip and tilt ofeither a point source or the barycentric shift of a compact diffusesource, the wide field Hartmann-Shack camera described herein attemptsto measure the geometrical distortion of an arbitrary extended source ofarbitrary dimensions. A point source and a compact diffuse source arespecial cases of the latter; therefore this wide field Hartmann Shackcamera is fully "downward compatible" to single point sources. Anotherdifference between the traditional H/S and this wide field H/S is thatthe former searches for a mono-axial measurement of tip-tilt centeredupon the barycenter of the beacon object, whereas the wide field H/Ssearches for the field variant tip-tilt across an 2 dimensional array offield angles. Given an undistorted image of an arbitrary beacon object,and a distorted image of the same, a methodology is outlined whereby thedistortion is measured and related to the tip-tilt of the localsub-aperture wavefront as a function of field angle. Yet anotherimportant advantage of the wide field wavefront camera, both the oneoutlined in this section and the ones outlined in the previous sectionsof this disclosure, is that they can function on zero-contrast objectsprovided that true anisoplanatism exists. Those practiced in the arthave loosely referred to this as the "blank wall" problem, or the "zerocontrast extended object" problem. Both the full field wavefront cameraof FIG. 7 and FIGS. 26A-B, and the wide field H/S camera of the presentsection, illustrate that this problem is solved due to the anisoplanaticsignal and can be measured by moving beyond the barycentric, monoaxialapproach to the Hartmann-Shack design.

Definite optical consideration not necessary in the traditional H/Ssensor come into play when performing detailed design of the wide fieldH/S camera. This disclosure points out and describes these extraconsiderations. After discussing these new considerations, the whollydifferent methodology of determining geometrical distortion will beoutlined and described in detail, then followed by a description of howthe measured geometrical distortion is transformed into a 2D array ofpupil function estimates.

New Optical Design Requirements

The general optical design of a wide field H/S camera is very similar tothe design of a traditional H/S sensor. In fact, it is easily like thatdepicted in FIG. 4 except that the field angle sampling aperture, 46, isreplaced by a simple square field stop that limits the field of view ofthe wide field H/S camera to some extent, such as 36 arc-seconds by 36arc-seconds. Furthermore, the determination of precise values of focallengths on the field lens, 50, the design variables on the lensletarray, 52, the bandpass characteristics on the spectral filter 38 anddichroic 28, the focal length on the magnificent optics 60, and the CCDpixel spacing specifications on the detector array, 10, may or will allbe modified by the extension to a wide field of view. Be this as it may,the elements, except for their detailed specifications as stated, remainthe same as in FIG. 4 and its corresponding written description. Ineffect, there will then be no true "field angle modulation" of theexiting light beam, 48, as there was in the previously disclosedwavefront camera designs. The measurement of geometric distortion itselfwill be the mechanism for distinguishing field angles and measuringtheir locate tip-tilt.

Measuring the Field Variant Tip-Tilt on and Extended Beacon

Reading from FIG. 4, each sub-aperture in the lenslet array, 52, will beprojecting low resolution image of the light distribution existing atthe field stop at the nominal focal plane, 66, onto a local region ofthe detector array, 10. This light distribution at the nominal focalplane, 66, is merely a broadband and generally chromatic and speckledistorted image of the beacon scene. The principle of the Hartmann-Shacksensor therefore still applies and thus a superposition of pupil images,weighted by the true brightness distribution of the celestial source, isprojected onto the lenslet array, 52, the superposition being one oflight from all the various field angles across the field stop. The imageprojected onto the detector array, 10, is therefore isolating the localtip-tilt of the incoming wavefront just as with a point source image,only now it is manifesting itself as a simple geometrical distortion ofan extended object due to the anisoplanatic changes in the localtip-tilt. This is to say, the wavefront tip-tilt at a given sub-aperturechanges as a function of field angle, and this change manifests itselfas a rather simple geometrical distortion on the detector array, ratherthan merely a barycentric shift of a point source.

It remains to determine as accurately as possible what the "undistortedbeacon image" looks like, akin to the need to find a nominal center orreference center in traditional H/S sensors, and once determining this,to then find the instantaneous distortion of a given exposure from this"reference image". Complicating the situation is the anisoplanaticbehavior of the amplitude of an incoming wavefront, i.e. thescintillation. We must also be explicit in now this additional effect ismeasured and dealt with. Once the geometric distortion between areference image and an instantaneous image is found, the localdifference value directly translates into local tip-tilt and thereforeall previous methods of traditional H/S sensors can be applied to obtaina full aperture pupil function for any given field angle which hassufficient beacon photons to produce a reasonable signal to noise ratio.

Determining the Reference Image for Each Sub-Aperture

The first step in the process is to determine a precise reference imageof the beacon object contained within the field stop at the nominalfocal plane 66 for each and every sub-aperture of lenslet 52. Thisreference object needs to be determined on the same pixel spacing scaleas the underlying detector array behind the lenslet array. A less purebut probably acceptable way is to determine a global reference imagewhich is then applied to each and every sub-aperture of the lensletarray 52. The former method of determining each separately merely canaccount for local inhomogeneities in not only the lenslet sub-aperturebut also the local detector array and other spatial non-uniformities inthe optics.

The reference image must have roughly equal, or slightly more, highspatial frequency signal level as any given distorted exposure. Thiswill help ensure higher signal to noise ratio measurement of thegeometric distortion. At the same time, this requirement complicates thederivation of the reference image. Essentially, the reference image isthe same image which would result if no distorting mechanism is present,i.e., if no atmosphere were present. There appear to be severalapproaches which can achieve the goal of optimally determining thereference image.

The preferred method of this disclosure is to design the wavefrontsensing optics in such a way that the individual MTF properties of agiven lenslet will be worse than the natural seeing conditions of theprimary telescope imaging system. In this way, a primary long exposureimage of the beacon object can be gathered on the primary detector andlow pass filtered to the same MTF of each sub-aperture and resampledonto the unique pixel grid behind each sub-aperture. It remains to takea long term exposure on the wavefront detector 10, and to align thecenters of gravity of the primary detector "reference image" with thelocal centers of gravity of the sub-aperture images. If there is a knowndegree of innate geometrical distortion existing in the lensletsthemselves, this systematic distortion should be added to the referenceimage as well.

Other methods of determining the reference image certainly may exist orbe developed in the future. Even using the long term sub-imagesthemselves is an alternative, but doing so will need some compensationwhen it comes time to measure the instantaneous geometrical distortinabout the reference image. MEASURING THE INSTANTANEOUS DISTORTIONBETWEEN THE REFERENCE IMAGE OF A GIVEN SUB-APERTURE AND AN ACQUIREDSHORT EXPOSURE SUB-APERTURE IMAGE ON DETECTOR 10

There are tow distinct approaches presented for measuring theinstantaneous distortion between the reference image and a given"distorted" image of a sub-aperture. No doubt there will be otherapproaches developed in the future, each with their own merits. It isimpossible to say which of those presented here is the preferredembodiment since each has its own merits relative to a cost/performancecomparison. The line integral method does argue to be the moreanalytical and possibly the more trustworthy, whereas the search methodhas conceptual simplicity and corresponding simplicity ofimplementation.

Description of the Problem

FIGS. 29A-D show an oversimplified illustration of an image and a purelydistorted, non-inverting version of that same image. The term "purelydistorted" refers to the idea that only a geometrical distortion, orwarping, has been applied to the undistorted image to produce thedistorted image. Non-inverting refers to the idea that if theundistorted image were made up of a large number of randomly directedsmall arrows, no single arrow would outright reverse direction in thedistorted image. FIG. 29 oversimplified in that the box would generallyhave a carved type distortion rather than a simple boxlike distortion,that is, the figure on the upper right would have curving lines.

The problem is to determine, as precisely as possible, the geometricaldistortion function (FIG. 29D) which maps the undistorted image (FIG.29A) to the distorted one (FIG. 29B). (FIG. 29C shows the distortionoverlaid on the original image.) Mathematically speaking, this is a twodimensional array of two dimensional distortion vectors. The problem tobe solved is really twofold: how to do this in simpler and ideal"noiseless" cases, and how to do this once noise, blurring,scintillation effects, and general error sources are added into thepicture.

Generic Pre-Processing

FIG. 30 contains an illustration of how the image of FIG. 29 might lookas seen by the detector 10 of FIG. 4. The scale is somewhat reduced andwe are assuming that there is only a 2 by 2 array of sub-aperturelenslets on the lenslet array 52. Hypothetical "distorted" images areshown on the left; the undistorted "reference image" is shown at theright. Note that the apparent geometrical distortion of one image to thenext changes. This is due to the anisoplanatic behavior of the localtip-tilt at each and every sub-aperture of the lenslet array 52 of FIG.4, and due to the difference of the tip-tilt between sub-apertures.Those practiced in the art will note that there are actually two primaryeffects taking place in the distorted images of FIG. 30: 1) thegeometric distortion of the image, and 2) a spatially varyingamplification and attenuation of the image due to anisoplanaticscintillation. The latter scintillation effect could not be easilydepicted in the graphics but should be understood to be present.

It turns out that measuring both of these distorting mechanismsindependently would be of use to the full distortion compensatingsystem. On the other hand, if all we wanted to recover was the puregeometric distortion, it would still be advantageous to separate the twoeffects.

It is clear that geometric distortion alone does not effect the totalbrightness of an object. Another way to state this is that the DC valueof the fourier transform is unaffected by geometric distortion. On theother hand, it will be found empirically that the scintillation functionfinds its maximum value, statistically, at the DC value of the fouriertransform and falls off rapidly as a function of spatial frequency. Thefact that these two effects occupy opposing regions of fourier spacewill therefore be exploited in the interests of separating the twoeffects.

Generally speaking, simple methods of separating the two can merelyassume that the two can be completely separated through a spatialfiltering process and that the resulting error can be simply tolerated.Moderately complicated methods can at least estimate and attempt tocompensate for the spatial-frequency cross-talk of the two effects,while ultimately, a full three dimensional model of the atmosphericphase distortion can be built which could form the basis for a VanCittert type iteration between ray-traced scintillation and phasedistortion (the derived function) and measured scintillation and phasedistortion (based directly on the data). Again and as always, we have avariety of sophistication we can apply to a given problem, each tradingoff cost with accuracy.

For present purposes we shall utilize a simple methodology, that ofmerely low-pass filtering both the reference image and each distortedimage, then adding the difference image between the two onto thedistorted image, producing a brightness-normalized distorted image withmuch of the scintillation effects removed (though some of the geometricdistortion effects will be removed in the process). As such, we findthat we have performed a gross separation of the two effects. Thedifference image, as a ratio to the reference image, is roughly ameasure of the field-variant scintillation function for each givensub-aperture, though both cross-talked geometric distortion andimprecise low-pass filtering will introduce errors. Thescintillation-normalized distorted image will be left largely with puregeometric distortion and will be input to the methodologies of the nextsections, through cross-talked scintillation and imprecise low-passfiltering will also introduce some errors.

Details on the exact steps involved in attempting to remove thescintillation function are contained in FIG. 31. Given a data image suchas that in FIG. 30, step 800 begins a loop through all sub-images of thereference object. Step 802 asks for a global registration between thecurrent distorted image and the reference image, and it suggests usingthe peak of the cross-correlation function between the two images.Indeed, inspection will show that for an image such as that depicted inFIG. 30, a cross-correlation might produce some strange results and notgive a clear peak in cross-correlation. It may then be prudent toperform a low pass filter on each image and repeat thecross-correlation. Step 804 then provides the options of removing knowngeometric distortion of the lenslet itself, if this refinement ispractical and worthwhile. Step 806 then performs a low pass filtering onboth reference image (it only needs to be performed for the firstsub-image then saved in a buffer, to be picky), and the distorted image.The suggested form of the filter is a gaussian, and the suggested FWHMis in the range of 5 to 10 arc-seconds, but these values ultimatelyshould be derived from empirical data which supports maximum separationof the two distortion signals. Step 808 then subtracts the filtereddistorted image from the filtered reference image, giving an estimate ofthe actual scintillation once it is normalized to the reference image.This estimate can then be used within the main processing routine ofFIG. 10 in the disclosure. Care must be taken, however, lest areas withno brightness whatsoever are encountered, wherein a weighting schemeshould be used based on the number of photons which are providing signalto the measurement of the scintillation. Step 810 acknowledges that thenormalized negative is a good estimate of the scintillation function.Step 812 then adds the difference image to the original (not filtered)distorted image to produce the "normalized" or "scintillation-removed"distorted image for the given sub-aperture of the main loop. Step 814closes the for-next loop begun in step 800.

Obtaining A Geometric Distortion Estimate Between A ReferenceUndistorted Image And A Distorted Image: 2 Approaches

The following two methods each produces an estimate of the geometricaldistortion obtaining between two generally identical, other thangeometric distortion, images. Having used the pre-processing methodabove of normalizing out the scintillation effects, we can now assumethat the error sources are simple noise and blurring and notscintillation. The problem here addressed has generally been describedin FIG. 29, where we will be searching for the bottom right diagram,given the upper left and the upper right as inputs.

The Iterative Fat-to-Thin Line Scan Method

FIG. 32 makes a basic point before diving into the main outline of thisgeneral approach. FIG. 32A is a simple one dimensional brightnessdistribution. FIG. 32B is the same distribution which has beengeometrically distorted. FIG. 32C shows the integrals of each, runningleft to right, overlaid on top of each other. The primary point of thisexercise is to note that the displacements of the ordinates between thetwo integral curves is a direct measure of the geometrical distortionbetween the two images--that is, in the case of pure non-invertingdistortion.

For one dimensional brightness distributions, therefore, it is a rathersimple operation to measure the distortion: merely plot out the graph asin FIG. 32C), and then directly generate the distortion plot (FIG. 32D)based on the abscissa differences of each and every ordinate, plotted onthe scale of the abscissas.

Moving on to the two dimensional case, the situation becomes morecomplicated. It is only mildly more complicated if the typicaldistortions are on the scale of one pixel or less, and rathercomplicated if the distortions are larger than one pixel and the imagehas high contrast, such as the case in FIGS. 29 and 30 (assuming thepixel spacings are much finer than the box spacings). Fortunately, inatmospheric distortion conditions the distortion is generally not sopronounced as depicted in FIGS. 29 and 30, and indeed is typically onthe order of one pixel in extent (depending on the specific opticaldesign, however). Thus, this first described method of the fat-to-thinprocess is worthy of consideration.

In a sense, the design of a high signal to noise ratio system preferslargely distortions because this will enable better precisionmeasurements of that distortion. This is to say that the desire forminimal distortion in the service of simplifying the method of findingthe correct distortion is in conflict with the desire for obtaining ahigher signal to noise ratio on finding that distortion.

When we apply a straightforward straight line integral to both images,the integral of the distribution along any given line on the distortedimage has brightness values from both sides of the line adding spurioussignal to the overall integration, thereby rendering the one-dimensionalmethod of FIG. 32 useless when applied uncritically. The supposedgeometric distortion being found can either be actual distortion oractual signal merely wrapped in from one side of the line or the other.

The line integral method, however, can be salvaged if we approach theestimation of the two dimensional distortion in an iterated fashion,first finding gross and very low spatial frequency distortion, thenrefining the spatial frequency estimation as we integrate along "curvedlines" which take into account the distortion profile on the orthogonalaxes and which have been roughly estimated from previous iterations.Ultimately this approach is one which must be justified on an empiricalbasis for any given turbulence situation.

FIG. 33 helps us take a pass through the required steps. The methodbegins (FIG. 33A) by taking broad line integrals along the orthogonalaxes, where broad refers to weighted spatial averaging along theperpendicular axis to the direction of the line integral path. Theresulting integral is further low-pass filtered and is compared to thesame process performed on the reference image. This process is appliedto a series of integrals first sweeping horizontal line integrals alongthe vertical axis, then sweeping vertical line integrals along thehorizontal axis. The line integral shown in FIG. 33A is of the formervariety. The total distortion found after a first full sweep is made isthen assumed to be a rough first guess at the distortion.

The spatial averaging weight function is then thinned up, and the newintegration path is taken along the estimated distortion lines from theprevious pass (FIG. 33B). In this way, a new distortion estimate isformed with slightly higher spatial frequency accuracy, while at thesame time following a better estimate of the proper integration path.After a few or several iterations, the method should converge upon thetrue geometric distortion (FIG. 33C). As with all iteration routines,convergence can be regulated by appropriate damping variables, such asupdating estimates by only one half their measured difference values,thus ensuing proper convergence.

(In FIG. 33, the thickness of the lines represent the width of theGaussian weighting functions along the orthogonal axis to the directionof travel of the line integral. The line integrals are first swept downthe y-axis as shown, then across the x-axis, together sweeping acrossthe entire two dimensional distortion. It may look simple shough, but wemust remember that there will be noise, error, and residualscintillation effects in real data, not to mention the potentialdifficulties of slow or local minimum convergence.)

The Global Search Method

The second method of estimating geometric distortion, the global searchmethod, is rather straightforward, but has the potential disadvantage oflarge computational loads and the difficulties associated with searchalgorithms including slow convergence times and the pitfalls of localminima. Nevertheless, search methods are also known to function properlyif these difficulties are well considered, and hence the followingtechnique is a viable one for determining the distortion function.

Generally, the unknown variables can be any set of functions whichadequately describe the two-axis shifting of an arbitrarily dense twodimensional array of control points. (i.e. reseau points). More simply,the unknowns can simply be the x-y offsets of some small grid of controlpoints, with inter-spacing points appropriately interpolated, such asthe bottom right drawing of FIG. 29.

The search can have as the primary metric the mean squared differencebetween the "redistorted" image, i.e. the distorted image geometricallywarped by the current search estimate vector, and the reference image.Pixels between the control points are, as mentioned, stretched viainterpolated amounts, and inter-pixel interpolation values are generallyuseful for more accurate results.

For the sake of speed, one pass of the "fat line integral" approach ofthe previous section can be used to get a grossly correct starting pointfor the search algorithm, and in the early stages only a certainpercentage of pixel values, say one in every four, need to be used tocalculate the mean squared error.

Though this search method is a bit more complicated than given creditfor in this summary, it is generally considered to be a mature prior artand the detailed implementation is a matter of typical refinement withinthat prior art.

Notes On An Anticipated Better Mode

There is a third method of determining the geometrical distortionbetween two images that is presently being worked on but has not yetbeen reduced to practice and verified. It is anticipated that thismethod will eventually be the preferred embodiment. It goes as follows.

The Line Integral Space Averaging Method

The method described here is a fully analytic method, to use this term abit loosely, in that it is a straightforward single iteration method offinding the precise geometric distortion function. Despite thisadvantage, it is strongly suggested that benchmark testing be performedon all three methods, or any newly devised ones, finding out whichmethod has the least computational overhead and ensuring that they allgenerate substantially similar estimates of the geometric distortion. Inorder to fully describe the line integral space method, some backgroundmust first be discussed.

FIG. 35a has an arbitrary line in two dimensional space. Provided thatthis line does not pass through the origin, it is uniquely identified bythe point which is closest to the origin. Lines passing through theorigin can be uniquely identified by the angle the make with the x axis.FIG. 35b formalizes this notion and creates an extended polar coordinatespace which uniquely defines all possible straight lines in a twodimensional euclidean space. At the origin, there is precisely 180degrees worth of unique lines, and therefore the r=0 point is expandedinto a circle, and the redundancy of line description in theta space isrepresented by a thick line on the upper half of the circle and a dashedline on the lower half.

FIG. 35c now defines a "linear integral transform" on this new extendedspace. For each point in the domain of the new extended space, the valueof a function is defined as the indefinite line integral of any giventwo dimensional function defined in normal 2D euclidean space, where theline integral is taken over the straight line specified by the point inthe domain. This is a unique mapping of a two dimensional function ontothe extended Line Integral Transform (LIT) space. In some sense, this isno different than any other functional transformation such as Fourier,LaPlace, or whatever.

(In more detail, FIG. 35C shows the definition of a Line IntegralTransform of a given two dimensional function f. F is a function of r,theta, which are the usual polar coordinates, plus the r=0, as definedabove. The integral is indefinite and is a line integral taken over thestraight line defined by r, theta. Generally speaking, spatially bound,continuous functions, f, will have a unique line integral transform. Theinverse transform is merely treating all points and values in F as aninfinite system of linear equations and solving for f. In spatiallydiscrete digital images, however, a finite set of equations willsuffice.)

Where is this taking us? Well, the first thing we need to realize isthat it is possible to transform back from the new space onto theconventional euclidean space. This can be done, to be overly pedantic,using a Hilbert-like transform where each point in the new space isexplicitly the addition of the function values along the line in theeuclidean space, thereby defining an infinite set of equations with anequally infinite set of unknowns, the classic uniquely mapped Hilbertformulation. Getting back to the simpler discrete world of digitalimages, it reduces to a large but nicely intuitive set of linearequations. For distortion problems in particular, where the distortionsare generally of low spatial frequency relative to the pixel spacing,this set of equations can become quite small and tractable.

We next define a new line integral called a radial line integral overthe transform space. The "radial" refers to the fact that this lineintegral has a constant theta in transform space, and only these typesof integrals will be utilized in the following steps. The integration isperformed on F, not f, and is performed in the Line Integral TransformSpace. Note that any given delta function in our original function firstbecomes a quasi-circle in LlT space (quasi because of the fat zero), andthat any given radial line integral will be a step function with thestep occurring at the appropriate point.

We then take a set of radial integrals equally spaced in theta space.Each integral is performed both on the reference image and on thedistorted image. The differences of the abscissas on any given ordinateis then interpreted as the integrated warp caused by the geometricaldistortion, just as it was in FIG. 32. For any given radial integral,these differences are found and plotted out along the radial lineitself. In this way, a new distribution is found plotted out in the lineintegral transform space. This is the line integral transform of thetrue geometrical distortion. It remains to set up this new distributionas a discrete set of equations and to solve this set of equations.

The resultant radial integration differences are then fed into a set oflinear equations in order to transform them back onto our normaleuclidean space. Solving for this matrix solves for the distortionfunction.

Some Final Notes On Inverting And Non-Inverting Distortion

It should now be clear that the preceding methods make use of theassumption of non-inverting distortion. This is not to say they will notfunction on inverting distortion, they will merely begin to exhibitdifferent and generally rapidly increasing errors in the results thatthey produce. It will be important for engineering considerations tobalance the needs of signal to noise ratios with the statistical onsetof these error sources deriving from inverting distortion, and to do soas a function of the prevailing strength of turbulence; which, asalready mentioned, is not much of a consideration for normal nighttimeatmospheric turbulence conditions.

Functioning On A Sparse Star Field Or An Image With Grossly Non-UniformContrast

A very important general case is now described, which also leads tofurther details on the implementation of the preceding methods. Thesalient features are drawn in FIG. 34. Essentially, the determination ofgeometric distortion is quite contingent upon either of both of A) localcontrast, or B) anisoplanatism. The two actually reinforce each other.That is, if both are present, the signal to noise ratio of themeasurement of the geometric distortion becomes that much better. Asimportantly, there must be sufficient detected photons in order to reachsome nominal signal to noise ratio performance.

Both methods described above must have modifications made whenconsidering objects which are either sparse or which have grosslynon-uniform contrast characteristics. Indeed, the modificationsdescribed herein should be at the core of the preferred embodiment whenusing any of these methods.

FIG. 34A shows a figure which has two randomly distributed point sourcesand two extended sources as an example. (The upper left extended objecthas lower contrast than the object at the lower right.) Inspection ofthe image shows that only at the places where there is an optical signalcan the measurement of geometric distortion take place. The implicationsof this are that only reseau points in the vicinity of the opticalsignal will get decent signal to noise ratios on their shiftmeasurements.

There are two basic ways of dealing with this: either use all the reseaupoints as if an entire S/N-ratio-homogenous image is present, or respacethe reseau points to "bunch around" the areas of maximum signal, so thatall reseau points have generally equal signal to noise ratios.Ultimately, both approaches will get to the same end point.

The point about the signal to noise ratio being contingent on either orboth of A) local contrast, and B) anisoplanatism is somewhat arcane butrelatively important. Conventional wisdom has it that zero-contrastextended objects do not give rise to a shift signal in wavefrontsensors, whereas FIG. 34, as well as the detailed discussion surroundingFIGS. 26A and 26B of the earlier disclosure, has shown thatanisoplanatism itself creates a shift signal even on zero-contrastobjects with at least some brightness. This fact alone is a performancedifferentiator between prior art mono-axial wavefront sensors and thewide field wavefront sensor of this section and the full field versionof a previous section.

Turning to a detailed look at FIG. 34B, the top line integral has awaveform of relatively constant positive slope, where it is labelled as`Anisoplanatic signal only`. This means that only the squeezing andstretching of the signal due to anisoplanatism will create a shiftsignal in this region. On the other hand, the lower line integral has`bumpy` signal in the extended object, so that both the squeezing andstretching of the signal due to anisoplanatism and the movement of theobject features will create a shift signal in this region. Thus, thesignal to noise ratio on the measurement of the geometrical distortionwill be higher on the high contrast extended object than on the lowcontrast extended object.

(FIG. 34B depict two line integrals showing how point sources and thehigher contrast extended object provide for high "shift signal" ofgeometric distortion, where the low contrast object still has "shiftsignal" but must rely purely on anisoplanatism to create the signal.)

Finally, the scintillation function measurement is almost exclusivelyrelated to just the primary local photon count and has little to do withthe local contrast. Thus, at the end of the day, the spatialnon-uniformities of the signal to noise ratios on the scintillationfunction and on the geometric distortion functions are only slightlycorrelated. This fact, however, is of little theoretical concern toultimate system performance.

Optimizing The Optical Design For The Wide Field Hartmann-ShackWavefront Camera

Generally speaking, it is important to vary the parameters of theoptical elements in order to maximize signal to noise ratios andminimize spatial filtering on pupil function estimates, all whileminimizing the number of photons required. FIG. 8 had a complete examplespecification for the full field wavefront sensor. Here we explore atypical example of specifications for the wide field Hartmann-Shackcamera.

The good news is that for most astronomical observations from nominalsites, inverting geometric distortion is quite rare. Empirical proof canbe readily obtained by merely imaging a binary star system and observingthat the anisoplanatism of the tilt--along the radial arm connecting thetwo images--is quite small and that the two stellar images effectivelyvibrate in unison. On very rare occasions, for our hypothetical site,the two images will merge or cross. The problem of inversion will becomea significant issue only when the ambient seeing conditions become quitepoor, such as r₀ =1 cm for example.

On to the more practical matters, a fundamental consideration in thedesign of the optical elements of the wide field Hartmann-Shack camerais trading off the signal to noise ratio on determining local geometricdistortion sensitivity. In general, the precise measurement of thegeometric distortion would like to see finer pixel spacings, while thedesire to collect photons into a single photosite for the purpose ofreaching dimmer beacons would dictate more coarse pixels "fine" and"coarse" refer to higher and lower magnifications respectively.

The "simple amplitude" of typical geometric distortion in r₀ =10 cmseeing is, in an rms sense, roughly in the one half arc-second region orless. Accurate data for any given target site, as a function of theseeing parameter, is easy to obtain by merely logging centroid positionsof stars and removing any tracking errors. Geometric distortion on theother hand can generally be measured to a decent inter-pixel level,which 1/10th of a pixel to 1/50th of a pixel is not untypical forreasonable light levels. Again this is in an rms sense.

Taking these very crude postulations, and positing once again the basictelescope arrangement of FIGS. 4 and 8, we can determine that a pixelspacing of 1 arc-second per pixel is reasonable for at least imagingbrighter objects such as planets. Here we find roughly a five totwenty-five signal to noise ratio on geometric distortion signal tonoise. On planets in particular, we must remember that we will havethousands of pixels giving us signal values, and the aggregate will thusgive us even better signal to noise performance due to a kind ofaveraging process. We can also target a 36 by 36 arc second square fieldof view for both the primary detector and the wide field Hartmann-Shackcamera, only because this will work out nicely for a 1024 by 1024 arrayon our 2.5 meter aperture telescope (Both Jupiter and Saturn can get abit bigger than this, and for those cases, the field of view would needto be enlarged using the basic principles herefollowing.

The basic optical elements are as they were in FIG. 4, except that theprimary array 8 should be increased to a 2048 by 2048 detector array inorder to image the larger field of view of 36" by 36", and the FieldAngle Sampling Sub-aperture Array, 46, is replaced by a Prime FocalPlane Field Stop, 47, which is merely a square aperture with dimensionsoutlined below. The incoming F/10 optical beam 32, FIG. 4, will notencounter and optional magnification stage, 60, as in FIG. 4, and thebeam will pass through the field lens 51, with a nominal focal length of134.3 mm, and will be occluded by a 4.5 mm by 4.5 mm square field stop47 (not shown, replaces Field Angle Sampling Sub-aperture Array, 46).The beam then fans out, passes through an open shutter, 42, during anexposure, and encounters the lenslet array 52 at the focal conjugate ofthe telescope pupil plane as imaged by field lens 51. The detectorarray, 10, is chosen as 1024 by 1024 in dimensions with a 15 micronsquare pixel spacing. The aperture size on each individual lenslet onthe lenslet array, 52, is chosen as 540 microns, with a individual focallength of 16.2 mm. The spacing from the first focal plane 66 to thelenslet array is a nominal (optical) 135 mm while the spacing from thelenslet array 52 to the detector 10 is a nominal 16.2 cm. Simplecalculations will show that this arrangement provides for an effectivesub-aperture coverage of the pupil function of 10 cm, or 25 lensletslinearly covering the 2.5 meter primary aperture. Moreover, it will befound that the plate scale at a given pixel on the detector array is at1 arc-second per pixel, or equivalently, one arc-second per 15 microns.These numbers are used as much for explanatory convenience as an attemptto optimize an actual design. Those practiced in the art will understandthat each design will determine its own unique specifications derivingfrom the constraints and trade-offs at hand. The sub-images of the fieldstop, 47, will be packed together with approximately four or five pixelsof border space between them.

A final note is in order concerning chromatic aberration of the lensletsystem. In the previous parts of this disclosure it was seen thatrelatively wide spectral band light beam is split into the wavefrontcamera. This is still the case with the currently described wide fieldHartmann-Shack (WFHS) camera. Now, however, the WFHS camera isexplicitly setting up geometrically distorted images of arbitrarilyextended objects rather than centroid shifts on point or highlylocalized sources. The chromatic aberrations will not change between thetwo systems, but the underlying processing methodologies may at somepoint being to experience a loss of signal amplitude due to thischromatic aberration. This effect should be noted in determining themaximum allowed spectral bandpass of the dichroic splitting arrangement28 and the optional filter 38, FIG. 4.

Converting The Distortion Measurements Into Field Angle Dependent PupilFunction Estimates

The preceding processing methods and their associated optical designlayout find both the scintillation function and the geometricaldistortion function of the sub-aperture images, functions themselves offield angle. It remains to convert these measurements into the formrequired by the processing methodology depicted in FIG. 10. What weactually require is to turn these measurements into the point spreadfunction estimates at the nominal field angle sample points depicted inFIG. 12A.

In general, this is exceedingly simple. The distortion function thusfound is a direct measure of the local tip-tilt, just as if a star waslocated precisely at each field angle sample point 306 of FIG. 12.Indeed, each field angle sample point has its local tip-tiltinterpolated for each and every sub-aperture of the lenslet 52, FIG. 4.Just as with the full field wavefront camera, each set of tip-tiltsacross the aperture is processed into an estimate of the pupil phasedistortion using the mature methodologies of the public domain. Thepossible and probable non-uniformity of the signal to noise level on thefield angle sample points is once again a fact of life and can merelyshow up as a spatial non-uniformity of the strehl ratios of final clearimages, 224, FIG. 10.

Now that we have a method for generating the field angle dependent pupilfunctions, we can replace the entire wide field Hartmann-Shack wavefrontcamera for the full field wavefront camera disclosed earlier and pick upwhere that disclosure left off. The full field camera may not beentirely obsolete, since some turbulence situations may be such that theoptics cannot reasonably produce good statistical non-inversion, and thefull field beam splitting design may offer some advantages in terms ofsimplifying and segregating the local potentially inverted distortions.This may also be of assistance in the quest for using dimmer and dimmerobject beacons, where "allowing" inversion may increase sensitivity. Bethis as it may, the simplicity and minimized optical surfaces of thewide field Hartmann-Shack arrangement should present a superior choicefor most applications.

Another Note On Anticipated Better Modes: Photon Counting Detectors AndCameras

The lower noise a detection system has, the better. The choice ofspecifying CCDs for both the primary detector 8 and the wavefrontdetector 10 is as much a matter of cost and simplicity as of ultimateperformance, barring cost considerations. Nevertheless, it is importantto be clear that the technology of spatial photon detection has been,and continues to be, evolving toward detection systems which some daywill be dominated by purely shot noise (photon arrival statistics,purely natural in deviation) and whereby instrumental noise is all buteliminated. The phrase and term `photon counting arry and camera` hasvariously been used to describe this ideal. It should be clear that theboth the detector 8 and the detector 10 will experience increasedperformance by using such a device as a photon counting detector. Anarrangement approaching this ideal is described in a copending U.S.application Ser. No. 08/171,562, filed Dec. 20, 1993, by Blouke andRhoads entitled "Wafer-Scale Photon Counting Array and Camera: AKA TheAvalanche CCD," incorporated herein by reference. (It should be notedthat the present application does not claim the technology disclosed inthe Blouke et al application; it is disclosed herein simply inaccordance with the inventor's duty in certain contracting states todisclose the contemplated best mode.). Other related arrangements areshown in Komobuchi, et al, "A Novel High-Gain Image Sensor Cell Based onSi p-n APD in Charge Storage Mode Operation," IEEE Trans. on ElectronDevices, 37:8, August, 1990, pp. 1861-1868, and Hynecek, "CCM--A NewLow-Noise Charge Carrier Multiplier Suitable for Detection of Charge inSmall Pixel CCD Image Sensors," IEEE Trans. on Electron Devices, 39:8,August 1992, 1972-1975.

It should be recognized that the particulars of the detector(s) employedin an embodiment of the present invention are entirely matters of designchoice for the system implementer; the present invention can beimplemented with a number of known detector technologies other thanthose described in the foregoing paragraph.

Derivative Versions of Resolution Enhancement Of Turbulence DegradedImages, For Use On Single Detector Telescopic Systems

So far, this disclosure has described an imaging system having twodetectors, one taking primary exposures and one taking wavefrontdistortion exposures. It so happens that a portion of the disclosedmethodologies can "trickle down" to systems which only employ a singleprimary detector array. In practical terms, the vast majority of imagingsystems have only a single detector (obviously), and thus a reasonableimprovement in resolving power on such systems, in the presence ofturbulence, might be of some use. In essence, the following method usesthe same exposure data both as primary data and as wavefront data. Byseparating out the two functions of the data, a methodology is presentedwhich should provide for reasonable resolution enhancement. Among manymitigating factors, it will be seen that the methods are best employedwhere an object scene does not change appreciably over, for example, aset of ten exposures. The reasons for this are that the methodspresented here work best when a multiple set of image data sets arepresent.

FIG. 36 gets us right to the heart of the methodology. FIG. 36A depictsa non-distorted image of an object; FIG. 36B depicts three examples ofthe same object with both geometrical distortion and spatially randomblurring. (Note that the blurring randomly roams around the image, bothacross spatial features and image to image.) The images on the right aretypical of, for example, short exposures (e.g. 10-100 milliseconds) on a30 cm to 1 m aperture telescope at a typical mountain site. It alsomight be typical of a terrestrial imaging system, of a somewhat smallerdiameter aperture and in somewhat higher turbulence conditions. Thepresent section of the disclosure describes a method whereby a set ofsuch images as on the right are taken and synthesized into a singleclear image of the object. The physical setup is therefore like that ofFIGS. 2 and 4, only now there is no dichroic splitter 28 and wavefrontcamera 11.

One novel aspect of the single detector arrangement disclosed herein isperhaps most clearly highlighted by comparison to certain prior artsystems. In particular, there are a wide variety of prior art imagingsystems which remove the global tip-tilt in a stream of video images.Specifically, there are cameras referred to as "steady-cams" which arewell known to remove "motion jitters." Likewise, cameras with eithermechanical and/or image processing technologies have been made which canproduce a steady image scene even when the camera is mounted on movingvehicles or aircraft. Finally, there are well known applications where aseries of "shifted images" are gathered and registered with each otherto effectively remove the "shift" between them, and then averagedtogether to form a single clear still image. The latter is typical ofplanetary imaging. In all of these prior arts, it is fair to say that"first order" global tip-tilt correction is being applied. A distinctionwith the arrangement described in this section is that the presentarrangement corrects beyond this first order.

In summary, the process described in this section performs a globalregistration of the distorted images, determines the local variations oftip-tilt in order to be removed (the geometric distortion), and thenpreferentially adds the areas of highest contrast via weightingcriteria. The weighting function is meant to keep track of localcontrast variations from one image to the next at a given location ofthe image, not from one part of the image to another. This acknowledgesthat often a given image can be crisp in one region and blurred inanother all in a single exposure, and that the crisp regions and blurredregions randonly wander around the image plane from one distorted imageto the next over a large series of E exposures. The point is that anygiven region will have its fair share of crisp data given enoughexposures.

In speckle imaging processing, "shift-and-add" techniques have been usedin certain applications, such as examining close double stars; a simplecentroid is determined, or a highest pixel value, the resulting x-yoffsets applied and the images added. In hindsight, the arrangement ofthe present section may be viewed as a rather involved extension on thiscore idea, applied to arbitrary extended objects. Likewise, it may seemfrom the literature that shifting and adding slightly offset images maybe used to reduce basic tracking errors and global tip-tilt in multipleexposure telescopic imaging. However, no specific description of such anapproach is presently known to the inventor, and this single detectorarrangement basically uses the global registration merely as a kind ofpre-processing of the images, or, to use analytic vernacular, the simpleremoving of the lowest order of geometric warp--global tip-tilt. It isgenerally acknowledged by those practiced in the art that these simpleshift-and-add techniques have their limitations and will only produce alimited resolution enhancement; thus the dual detector set-up of theearlier disclosure is believed necessary to achieve near diffractionlimited performance on larger aperture telescopes and in more turbulentconditions.

FIG. 37 steps through the single detector enhancement methodology. Webegin once again by gathering a series of short exposure digital imagedata sets of the object, step 900. The next step, 902, is to calibratethese digital images using the very common techniques of offset and gainnormalization to obtain linear measurement values. This is also known asdark framing and flat fielding. The next step 904 is to register allimages onto one global x-y offset, effectively removing tracking errorand global tip-tilt of the images. This process is well known in the artof image processing and consequently has several approaches. One suchapproach is to cross-correlate images and find where therecross-correlation peaks. Certain provisions for rotation of images mayneed addressing in systems which do not have good mechanical trackingand rotating properties.

Once two given images are "registered", it remains to resample one imageonto the grid spacing of the other. This too is a familiar imageprocessing technique and reference is made to standard texts dealingwith these basics. The images can then be effectively averaged together.The registered images now become the working exposures for the nextsteps.

Step 906 averages all E exposures together pixel by pixel to create the"standard image" which we explicitly acknowledge is a slightly lowerresolution image than what we are ultimately after: we then apply a lowpass filter to the average image via step 908. The filter is typicallyof gaussian shape with a full width half maximum of 5 arc-seconds or so,but in any event is subject to experimentation. This resulting image iscalled the standard scintillation image. Step 910 then steps througheach and every exposure, filters it (step 912) with the same filterapplied to get the standard scintillation image, and adds the differencebetween the standard and the newly filtered version to the originalexposure version, giving each exposure a "scintillation-correct version"(via steps 914 and 916), i.e. some large statistical amount of theanisoplanatic scintillation has now been removed from each exposure.

Step 918 now calls the earlier-detailed procedure which finds thegeometric distortion between the "scintillation-corrected image" and the"standard image" of step 906. Of note is that the standard image shouldbe generally of a lower frequency content than the scintillationcorrected images, i.e. it is a bit more blurred on average, and thatduring the geometric distortion mensuration process step 918, theresulting geometric distortion should contain a very high frequencyspurious signal that should be filtered out as in step 920. Finally,step 922 creates a new "undistorted image" by applying the reversegeometric distortion to the scintillation corrected image of step 916.

Step 924 next creates a "weighting image" which essentially is a measureof the local contrast, or "detail," of the geometrically corrected imagefrom step 922. The specifics of generating this image are quiteflexible. The preferred embodiment creates this image by, at each andevery pixel, assigning the new value of the local absolute value of theaveraged second derivative at that pixel location. This is generallycomputed on a floating point basis, later to be scaled to, for example,an 8-bit scale where the value of zero is mapped to the digital value of1 and the highest anticipated second derivative value is scaled to thedigital number 200. If the highest value exceeds 200 or even 255, thescaling must be changed; likewise if the highest value is well under thevalue 200. This image is then low pass filtered with generally the sameFWHM of the gaussian used to create the scintillation-corrected images.Once this weight image is finished, step 926 now steps through pixel bypixel then adding its "correct image value" weighted by its "weightingimage value" into an accumulating weighted average for each and everypixel of the final clear image. Step 928 completes the loop 910 suchthat when all E exposures are cycled through, step 930 then begins theloop through all the weighted average numerators and denominators ofeach and every pixel, performs the required division in step 932, closesthe loop 934, and out pops an optimally enhanced image estimate 936. Aswith any image processing routine such as the one outlined here, certainresidual effects may become apparent, and it is always best to specify apost-filtering process step 938 of a generic kind, with the hopes ofsmoothing out any undue anomalies.

A few more ancillary notes are in order. First, short exposures aregenerally advised so long as they don't trade off other performancegoals such as overall noise. The short exposures tend to "freeze" bothatmospheric distortion and minor tracking error. If either of theseeffects last too long and are integrated into one given exposure, theytend to irreversibly destroy high frequency information in the rawexposures. Second, narrow band filtering can also be effectivelyemployed. The filtering merely minimizes chromatic blurring onmonochromatic detectors, again assisting in retaining high frequencyinformation in the raw exposures.

Conclusion

This disclosure has aptly focussed upon an "optical telescope" operatingin the open atmosphere as the preferred embodiment of the invention, andhas numerous references to alternative embodiments. "Optical" has bothexplicitly and implicitly been used in its more general sense asreferring to the entire electromagnetic spectrum. The wavelength of 500nanometers was repeatedly used as an example central wavelength aboutwhich a narrow band of light is imaged. To use this wavelength as anexample one final time, it can be seen that whereas prior art systemssuch as adaptive optics systems and speckle imaging systems could indeedachieve near diffraction limited performance at this wavelength andothers like it, they could do so only within a steradian field coverageof 5 arc seconds in atmospheric conditions in the r₀ =10 cm. Thisdisclosure, on the other hand, specified how to attain near diffractionlimited imagery to arbitrarily wide fields of view (where again, an 8 by8 field angle sampling was merely an example). Furthermore, prior arttelescopes not utilizing either prior art speckle methods or adaptiveoptics methods could generally not resolve objects to better than 1.0arc seconds over any field of view (at 500 nanometers).

It should once again be noted that, while the foregoing disclosure hasaddressed in considerable detail certain telescope applications, theteachings thereof are more generally applicable to any imaging systemdegraded by temporally or spatially varying distortion.

The issue of alternative embodiments to those detailed herein has beenquite a recurring theme in this disclosure. Those skilled in the art oftelescope or imaging system design will be able to clearly recognize therange of design options available for creating other embodiments. Thoughthis principle is almost a cliche in the spirit of patent protection, itis nevertheless quite important to recognize this freshly for inventionswhich must perforce borrow upon not only decades but centuries of priorart; inventions which in some embodiments are best described as verylarge systems. It was with this in mind that FIGS. 1, 5, and 11A weredrawn and described with certain novelties highlighted: And it was withthis mind that so many design options were positively "encouraged"rather than merely being obliquely referenced. Therefore, I claim as myinvention all such embodiments as may come within the scope and spiritof the following claims and equivalents thereto.

What is claimed is:
 1. Apparatus comprising:a detection system forsensing wide field atmospheric distortion; and a processing system,coupled to the detection system, for generating wide field imagingcompensation data corresponding to said distortion.
 2. Acomputer-readable storage medium, having instructions stored thereon forcausing a computer programmed thereby to perform the followingmethod:reading a set of distorted image data; characterizing a fieldvariant transfer function; and adding to the distorted image data anensemble of locally corrected basis functions to compensate for saiddistortion, said ensemble corresponding to the characterized fieldvariant transfer function.
 3. The invention of claim 2 wherein thecharacterizing step includes searching on corrected image data toestimate the field variant transfer function.
 4. The invention of claim2 in which the instructions cause computer programmed thereby to:breakdown a set of corrected image data into an overdetermined set of saidbasis functions; perform locally applicable distortion reversion methodson local sub-groups of said overdetermined set; and accumulatedistortion reversed local sub-groups into a clear image.
 5. Theinvention of claim 2 in which the instructions cause the computerprogrammed thereby to characterize a plurality of field variant transferfunctions, each for a different set of distorted image data, and forminga weighted average of said distorted image data sets using signal tonoise preserving basis functions in accordance with said characterizedsets of field variant transfer functions to produce a corrected imagedata set.
 6. A computer-readable storage medium, having instructionsstored thereon for causing a computer programmed thereby to perform thefollowing method:receiving a set of detector data corresponding to atwo-dimensional representation of an extended source as imaged through adistorting medium; processing said set of detector data to discerngeometric distortion associated with the distorting medium; anddetermining, from said discerned geometric distortion, athree-dimensional turbulence structure of the distorting medium.
 7. Theinvention of claim 6 in which the instructions cause the computerprogrammed thereby to:receive a plurality of sets of said detector data;and processing said plurality of sets to discern said geometricdistortion.